From Risk to Resilience: Ecological Exploration of Psychopathology and Flourishing in Neurodivergent Children Facing Adversity

Project Overview

---

This study investigates the relationship between adverse childhood experiences (ACEs) and mental health outcomes in children with neurodevelopmental disorders (NDDs), such as Autism Spectrum Disorder (ASD) and ADHD. The study specifically examines the moderating role of child flourishing behaviors—curiosity in learning, emotional control, and task persistence—on the effects of ACEs on mental health outcomes.

Study Design:
The study utilizes data from the National Survey of Children’s Health (NSCH), focusing on children aged 6-17 with reported neurodevelopmental issues (N = 54,528, M_Age = 12.1). This cross-sectional analysis assesses the exposure to ACEs and its impact on three key mental health outcomes: anxiety, depression, and behavioral problems. Interaction terms in logistic regression models are used to assess whether child flourishing behaviors buffer the effects of ACEs on these outcomes.

Data Collection:
Key variables collected include:

• Adverse Childhood Experiences (ACEs): Derived from parent-reported experiences such as exposure to violence, household dysfunction, and discrimination.
• Mental Health Outcomes: Parent-reported diagnoses of anxiety, depression, and behavioral problems, confirmed by health care providers or educators.
• Child Flourishing Metrics: Based on parent-reported behaviors such as curiosity in learning, emotional control, and task persistence, which were used to develop a Child Flourishing Index (CFI).

Analytic Approach:
Logistic regression models were employed to evaluate the association between ACEs and mental health outcomes, including interaction terms to investigate the moderating role of child flourishing. The models incorporated a comprehensive set of covariates, including child age, sex, race/ethnicity, socioeconomic status, caregiver mental health, neighborhood conditions, and healthcare access metrics.

Goals:
The study aims to:

• Examine the dose-response relationship between ACE exposure and mental health outcomes among neurodivergent children.
• Assess the protective role of child flourishing behaviors, with the goal of identifying resilience-promoting factors that can inform interventions for this vulnerable population.

Tip

For detailed information on the NSCH-derived variables, refer to the NSCH Data Dictionary.

For a primer on Neurodiversity: What is Neurodiversity?

Data Frame Initialization

---

Set up the R environment by configuring the CRAN repository, installing essential packages, setting the base path, and loading the primary dataset into the local environment.

# Set CRAN repository for consistent package installation
options(repos = c(CRAN = "http://cran.rstudio.com/"))

# Load or install the pacman package for efficient package management
if (!require(pacman)) install.packages("pacman")

# Use p_load to load or install the required packages
p_load(
  dplyr,           # Data manipulation
  tidyr,           # Data reshaping (e.g., pivot_longer)
  readr,           # Reading data (CSV files)
  knitr,           # Creating tables for reports
  kableExtra,      # Enhancing `knitr::kable()` tables
  ggplot2,         # Data visualization
  sjlabelled,      # Labelled data utility functions
  sjmisc,          # Additional data manipulation functions
  survey,          # Complex survey analysis
  lmtest,          # Testing models (Breusch-Pagan test) and robust standard errors
  sandwich,        # Calculating robust standard errors
  gtsummary,       # Creating regression summary tables
  ggeffects,       # Creating prediction plots
  gt,              # Adjusting frequency/contingency tables
  psych,           # Factor analysis
  crosstable,      # Cross-Tabulations
  ggpubr,          # Summary statistics + visualizations
  tidyverse,       # Data visualization supplement
  ggdist           # Distribution visualizations
  )

# Specify the 'base_path' where you can find your data files, ASSUMING they're in the same directory, & set it as WD
base_path <- "~/GitHub/Adversity-Neurodiversity/data_files"
setwd(base_path)

# Load primary data file
neurodivergent_data <- read_csv("~/Downloads/nsch_2023e_topical_SAS/neurodivergent_data_2016to2023.csv")

1

See details under the Data Extraction & Filtering (Archived Code) callout below!

Data Extraction & Filtering (Archived Code)

The file being used for these analyses is a subset of a “master” file, Data2_2016to2022.csv, containing the compiled contents of the data releases from the National Survey on Child Health’s 2016-2022 cycles. The master file is >500 MB (0.5 GB) in size, making it cumbersome to process in real-time. For efficiency and ease of analysis, I am using a smaller subset that contains only the relevant data. This subset was generated using the data extraction and filtering process detailed in the archived code below, in an effort to maintain transparency and ensure reproducibility of the workflow.

Additionally, data from the NSCH 2023 was also later merged into the subset file being used for this analyses. The code for the extraction/wrangling process is shown below.

# Define the path to the dataset
data_path <- "~/Downloads/Data2_2016to2022.csv"  
Data2_2016to2022 <- read.csv(data_path)

# Select specific variables, create 'neurodivergent' variable, & filter data based on age and neurodivergent status
neurodivergent_data <- Data2_2016to2022 %>%
  select(year, fpl, fpl_mean, sc_age_years, sc_hispanic_r, sc_race_r, sc_sex, higrade_tvis,
         k2q35a, k2q35c, k2q38a, k2q38c, k2q30a, k2q30c, k2q36a, k2q36c, k2q60a, k2q60c, 
         k2q37a, k2q37c, downsyn, downsyn_desc, k2q31a, k2q31c, k2q61a, cerpals_desc, 
         k2q33a, k2q33b, k2q33c, k2q32a, k2q32b, k2q32c, k2q34a, k2q34b, k2q34c, 
         ace1, ace3, ace4, ace5, ace6, ace7, ace8, ace9, ace10, ace11, ace12, k6q71_r, 
         k7q85_r, k7q84_r, talkabout, wktosolve, strengths, hopeful, k8q21, k8q30, fwc, 
         k2q42a, k2q42c, a1_menthealth, a2_menthealth, k2q01, k4q27, k4q28x01, k4q28x04, 
         appointment, available, issuecost, notelig, notopen, transportcc, family, family_r,
         k10q11, k10q12, k10q13, k10q14, k10q20, k10q22, k10q23, k10q40_r, k10q41_r,
         k10q30, k10q31, k9q96, sc_cshcn, k3q22, k3q20, instype, treatneed, k4q22_r,
         menbevcov) %>%
  # Creating 'neurodivergent' variable as NSCH does not explicitly inquire about neurodivergent status.
  # Using a list of diagnoses provided by NSCH to define neurodivergent individuals.
  mutate(neurodivergent = if_else(k2q35a == 1 | k2q38a == 1 | k2q36a == 1 | k2q60a == 1 | 
                                  k2q37a == 1 | k2q30a == 1 | downsyn == 1 | 
                                  k2q31a == 1 | k2q61a == 1 | k2q42a == 1, 1, 0)) %>%
  # Filter data to include only individuals aged 6 or above and identified as neurodivergent
  filter(sc_age_years >= 6, neurodivergent == 1)

# Resulting filtered data to be used for further analysis
write.csv(neurodivergent_data, "neurodivergent_data_2016to2022.csv")

# Define the path to the dataset
data_path <- "~/Downloads/nsch_2023e_topical_SAS/nsch_2023e_topical.csv"  
nsch_data <- read.csv(data_path)

# Create a new variable 'fpl_mean' representing the average of specific columns
nsch_data <- nsch_data %>%
  mutate(fpl_mean = rowMeans(select(., fpl_i1, fpl_i2, fpl_i3, fpl_i4, fpl_i5, fpl_i6), na.rm = TRUE))

# Select specific variables, create 'neurodivergent' variable, & filter data based on age and neurodivergent status (manually added in blank columns fpl and downsyn_desc)
# `menbevcov` was not collected for 2023
neurodivergent_data_23 <- nsch_data %>%
  select(year, fpl, fpl_mean, sc_age_years, sc_hispanic_r, sc_race_r, sc_sex, higrade_tvis, 
         k2q35a, k2q35c, k2q38a, k2q38c, k2q30a, k2q30c, k2q36a, k2q36c, k2q60a, k2q60c, 
         k2q37a, k2q37c, downsyn, downsyn_desc, k2q31a, k2q31c, k2q61a, cerpals_desc, 
         k2q33a, k2q33b, k2q33c, k2q32a, k2q32b, k2q32c, k2q34a, k2q34b, k2q34c, 
         ace1, ace3, ace4, ace5, ace6, ace7, ace8, ace9, ace10, ace11, ace12, k6q71_r, 
         k7q85_r, k7q84_r, talkabout, wktosolve, strengths, hopeful, k8q21, k8q30, fwc, 
         k2q42a, k2q42c, a1_menthealth, a2_menthealth, k2q01, k4q27, k4q28x01, k4q28x04, 
         appointment, available, issuecost, notelig, notopen, transportcc, family_r,
         k10q11, k10q12, k10q13, k10q14, k10q20, k10q22, k10q23, k10q40_r, k10q41_r,
         k10q30, k10q31, k9q96, sc_cshcn, k3q22, k3q20, instype, treatneed, k4q22_r) %>%
  # Creating 'neurodivergent' variable as NSCH does not explicitly inquire about neurodivergent status.
  # Using a list of diagnoses provided by NSCH to define neurodivergent individuals.
  mutate(neurodivergent = if_else(k2q35a == 1 | k2q38a == 1 | k2q36a == 1 | k2q60a == 1 | 
                                  k2q37a == 1 | k2q30a == 1 | downsyn == 1 | 
                                  k2q31a == 1 | k2q61a == 1 | k2q42a == 1, 1, 0)) %>%
  # Filter data to include only individuals aged 6 or above and identified as neurodivergent
  filter(sc_age_years >= 6, neurodivergent == 1)

# Resulting filtered data to be used for further analysis
write.csv(neurodivergent_data_23, "neurodivergent_data_2023.csv")

## Manually removed index column in both 2016-2022 and 2023 files 
## Manually add empty `family` & `menbevcov` columns to 2023 file
# Read in update 2023 data file
neurodivergent_data_23 <- read.csv("~/Downloads/nsch_2023e_topical_SAS/neurodivergent_data_2023.csv")

# Read in prior years of neurodivergent data
neurodivergent_data_16to22 <- read.csv("~/Downloads/nsch_2023e_topical_SAS/neurodivergent_data_2016to2022.csv")

# Merge two data frames with identical columns
neurodivergent_data_16to23 <- rbind(neurodivergent_data_16to22, neurodivergent_data_23)

# Write the merged data frame into a CSV file
write.csv(neurodivergent_data_16to23, "neurodivergent_data_2016to2023.csv", row.names = FALSE)

Analytic Data Preparation

---

Child Demographic Variables

Recode survey year into a factor, recode child sex as a categorical variable, recode race and ethnicity into descriptive categories, and create a combined race-ethnicity variable.

# Recode survey year into a factor variable
neurodivergent_data <- neurodivergent_data %>%
  mutate(survey_year = factor(year))

# Recode sex and race/ethnicity categories
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Recode sex into categorical factors
    sc_sex_cat = factor(sc_sex, levels = c(1, 2), labels = c("Male", "Female")),
    
    # Recode race categories into descriptive labels
    race = case_when(
      sc_race_r == 1 ~ "White",
      sc_race_r == 2 ~ "Black",
      sc_race_r == 3 ~ "American Indian or Alaska Native",
      sc_race_r %in% 4:5 ~ "Asian, Native Hawaiian, or Pacific Islander",
      sc_race_r %in% 6:7 ~ "Other or Mixed Race",
      TRUE ~ NA_character_
    ),
    
    # Recode ethnicity
    ethnicity = case_when(
      sc_hispanic_r == 1 ~ "Hispanic/Latino",
      sc_hispanic_r == 2 ~ "Non-Hispanic/Latino",
      TRUE ~ NA_character_
    ),
    
    # Create combined race and ethnicity category
    sc_race_cat = case_when(
      race == "White" & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic White",
      race == "Black" & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic Black or African American",
      race == "American Indian or Alaska Native" & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic American Indian or Alaska Native",
      race == "Asian, Native Hawaiian, or Pacific Islander" & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic Asian, Native Hawaiian, or Pacific Islander",
      race == "Other or Mixed Race" & ethnicity == "Non-Hispanic/Latino" ~ "Non-Hispanic Other or Mixed Race",
      ethnicity == "Hispanic/Latino" ~ "Hispanic/Latino of Any Race",
      TRUE ~ NA_character_
    ),
    sc_race_cat = factor(sc_race_cat, levels = c(
      "Non-Hispanic White", 
      "Non-Hispanic Black or African American", 
      "Non-Hispanic American Indian or Alaska Native", 
      "Non-Hispanic Asian, Native Hawaiian, or Pacific Islander", 
      "Non-Hispanic Other or Mixed Race", 
      "Hispanic/Latino of Any Race"
    ))
  )

Healthcare-Related Variables

Create healthcare-related indices including Healthcare Access, Barriers to Access, Insurance Type, and Insurance Adequacy. Recode variables for consistent formatting, compute indices, and categorize healthcare access levels based on defined thresholds.

# Create a Healthcare Access Index and categorize it
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    healthcare_access_index = case_when(
      k4q27 == 2 & k4q22_r == 1 ~ 2,  # High access
      k4q27 == 2 & k4q22_r == 3 ~ 1,  # Medium access
      k4q27 == 1 | k4q22_r == 2 ~ 0,  # Low access
      TRUE ~ NA_real_ # Assign NA for missing or undefined values
    ),
    HAI_category = factor(healthcare_access_index, 
                          levels = c(2, 1, 0), 
                          labels = c("High Access", "Medium Access", "Low Access"))
  )

# Recode 'treatneed' to standardize across survey years
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    treatneed_score = case_when(
      survey_year %in% 2018:2023 & treatneed == 1 ~ 0,  # Not difficult (2018-2023)
      survey_year %in% 2018:2023 & treatneed == 2 ~ 1,  # Somewhat difficult (2018-2023)
      survey_year %in% 2018:2023 & treatneed == 3 ~ 2,  # Very difficult (2018-2023)
      survey_year %in% 2018:2023 & treatneed == 4 ~ 3,  # Not possible (2018-2023)
      survey_year %in% 2016:2017 & treatneed == 1 ~ 0,  # Not a problem (2016-2017)
      survey_year %in% 2016:2017 & treatneed == 2 ~ 1,  # Small problem (2016-2017)
      survey_year %in% 2016:2017 & treatneed == 3 ~ 2,  # Big problem (2016-2017)
      TRUE ~ NA_real_  # Assign NA for missing or undefined values
    )
  )

# Create the Barriers to Access Index
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    barriers_score = rowSums(select(., appointment, available, issuecost, notelig, notopen, transportcc, treatneed_score) == 1, na.rm = TRUE),
    barriers_cat = cut(barriers_score, 
                       breaks = c(-1, 0, 2, 4, 7), 
                       labels = c("No Barriers", "Few Barriers", "Moderate Barriers", 
                                  "Many Barriers"))
  )

# Include insurance type as a categorical variable
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    insurance_type = case_when(
      instype == 1 ~ "Public Only",
      instype == 2 ~ "Private Only",
      instype == 3 ~ "Private + Public",
      instype == 5 ~ "Uninsured",
      instype == 4 ~ NA_character_  # Excluded 'insurance type unspecified', n = 85
    ),
    insurance_type = factor(insurance_type, 
                            levels = c("Private Only", "Private + Public", 
                                       "Public Only", "Uninsured"))
    )

# Create an Insurance Adequacy Index
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Assign scores to the categorical responses
    k3q20_score = recode(k3q20, `1` = 4, `2` = 3, `3` = 2, `4` = 1), 
    k3q22_score = recode(k3q22, `1` = 4, `2` = 3, `3` = 2, `4` = 1),
    menbevcov_score = recode(menbevcov, `1` = 4, `2` = 3, `3` = 2, `4` = 1, `5` = NA_real_), 
    # Excluded `menbevcov` = 5 'child does not need/use M/BH services', n = 15,310
    
    # Sum scores to create an Insurance Adequacy Index
    insurance_adequacy_index = rowSums(across(c(k3q20_score, k3q22_score, menbevcov_score)), na.rm = TRUE),

    # Categorize the Insurance Adequacy Index
    IAI_category = case_when(
      insurance_adequacy_index >= 10 & insurance_adequacy_index <= 12 ~ "High Adequacy",
      insurance_adequacy_index >= 7 & insurance_adequacy_index < 10 ~ "Moderate Adequacy",
      insurance_adequacy_index >= 3 & insurance_adequacy_index < 7 ~ "Low Adequacy",
      TRUE ~ NA_character_  # Assign NA for missing or undefined values
    ),
    IAI_category = factor(IAI_category, levels = c("High Adequacy", "Moderate Adequacy", "Low Adequacy"))
  )

1

k4q27 = Needed Health Care Not Received; k4q22_r = Mental Health Professional Treatment

2

treatneed = Mental Health Professional Treatment - Access Difficulty; appointment = Needed Health Care Not Received Due to - Getting Appointment; available = Needed Health Care Not Received Due to - Not Available; issuecost = Needed Health Care Not Received Due to - Cost; notelig = Needed Health Care Not Received Due to - Not Eligible; notopen = Needed Health Care Not Received Due to - Office Not Open; transportcc = Needed Health Care Not Received Due to - Getting Transportation

3

k3q20 = Benefits Cover Services; k3q22 = Allow to See Provider; menbevcov = Cover Mental Behavioral Needs

Family & Caregiver Context Variables

Recode caregiver education, family income, and caregiver composition variables to categorical values. Create a composite score for caregiver mental health, then categorize it into meaningful health levels.

# Recode caregiver education and family income categories
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Recode caregiver education into categories
    caregiver_education = case_when(
      higrade_tvis == 1 ~ "Less Than High School",
      higrade_tvis == 2 ~ "High School/Vocational",
      higrade_tvis == 3 ~ "Some College/Associate Degree",
      higrade_tvis == 4 ~ "College Degree/Higher",
      TRUE ~ NA_character_ # Assign NA for missing or undefined values
    ),
    caregiver_education = factor(caregiver_education, levels = c("Less Than High School", "High School/Vocational", "Some College/Associate Degree", "College Degree/Higher")),

    # Recode family income into federal poverty level categories
    fpl_category = case_when(
      fpl %in% c(50:99) | fpl_mean < 100 ~ "<100%",
      fpl %in% c(100:199) | fpl_mean < 200 ~ "100%-199%",
      fpl %in% c(200:299) | fpl_mean < 300 ~ "200%-299%",
      fpl %in% c(300:399) | fpl_mean < 400 ~ "300%-399%",
      fpl %in% c(400:999) | fpl_mean < 999 ~ "≥400%",
      TRUE ~ NA_character_ # Assign NA for missing or undefined values
    ),
    fpl_category = factor(fpl_category, levels = c("<100%", "100%-199%", "200%-299%", "300%-399%", "≥400%"))
  )

# Simplify caregiver composition based on family structure
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    caregiver_composition = case_when(
      family %in% c(1, 2, 3, 4) | family_r %in% c(1, 2, 3, 4) ~ "Dual-caregiver household",
      family %in% c(5, 6, 7, 8, 9) | family_r %in% c(5, 6, 7, 8) ~ "Single-caregiver household",
      TRUE ~ NA_character_ # Assign NA for missing or undefined values
    ),
    caregiver_composition = factor(caregiver_composition, levels = c("Dual-caregiver household", "Single-caregiver household"))
  )

# Create a composite variable for caregiver mental health
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    caregiver_mental_health_score = rowMeans(select(., a1_menthealth, a2_menthealth), na.rm = TRUE),
    caregiver_mental_health = cut(
      caregiver_mental_health_score,
      breaks = c(1, 1.5, 2.5, 3.5, 4.5, 5),
      labels = c("Excellent", "Very Good", "Good", "Fair", "Poor"),
      include.lowest = TRUE
    )
  )

1

Scores of caregiver mental health are averaged to account for single respondents (i.e., single-caregiver households).

Neighborhood Context Variables

Recode neighborhood binary items to 0/1 format and standardize ordinal neighborhood variables by subtracting their mean and dividing by the standard deviation.

# Recode neighborhood condition-related binary items to 0/1
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    across(c(k10q11, k10q12, k10q13, k10q14, k10q20, k10q22, k10q23, k9q96), ~ ifelse(. == 1, 1, 0))
  )

# Standardize neighborhood condition-related ordinal variables
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    across(c(k10q40_r, k10q41_r, k10q30, k10q31), ~ scale(., center = TRUE, scale = TRUE))
  )

1

k10q11 = Neighborhood - Sidewalks or Walking Paths; k10q12 = Neighborhood - Park or Playground; k10q13 = Neighborhood - Recreation Center; k10q14 = Neighborhood - Library or Bookmobile; k10q20 = Neighborhood - Litter or Garbage; k10q22 = Neighborhood - Poorly Kept or Rundown Housing; k10q23 = Neighborhood - Vandalism; k9q96 = Other Adult Child Can Rely On For Advice

2

k10q40_r = Child is Safe In Neighborhood; k10q41_r = Child Is Safe at School; k10q30 = People In Neighborhood Help Each Other Out; k10q31 = Watch Out for Other’s Children

Factor Analysis - Neighborhood Context Scores

Conduct factor analysis on neighborhood context variables to derive factor scores for physical resources, environmental threats, and collective efficacy, and add these scores to the primary data frame. Remove intermediate data frames from the global environment.

The factor structure was adapted from factor conceptualization detailed in Fan & Chen (2012), “Family functioning as a mediator between neighborhood conditions and children’s health: Evidence from a national survey in the United States.” Principal Axis Factoring (PAF) was used as the factoring method, and factor scores were calculated using the regression method, thus consistent with the approach in the referenced study.

# Select the items for each neighborhood factor
physical_resources_items <- neurodivergent_data %>%
  select(k10q11, k10q12, k10q13, k10q14)

environmental_threats_items <- neurodivergent_data %>%
  select(k10q20, k10q22, k10q23, k10q40_r, k10q41_r)

collective_efficacy_items <- neurodivergent_data %>%
  select(k10q40_r, k10q41_r, k10q30, k10q31, k9q96)

# Conduct factor analysis for each factor using Principal Axis Factoring (PAF)
fa_physical_resources <- fa(physical_resources_items, nfactors = 1, fm = "pa")
fa_environmental_threats <- fa(environmental_threats_items, nfactors = 1, fm = "pa")
fa_collective_efficacy <- fa(collective_efficacy_items, nfactors = 1, fm = "pa")

# Calculate factor scores for physical resources using regression method
physical_resources_scores <- factor.scores(physical_resources_items, fa_physical_resources, method = "regression")$scores[, 1]

# Calculate factor scores for environmental threats using regression method
environmental_threats_scores <- factor.scores(environmental_threats_items, fa_environmental_threats, method = "regression")$scores[, 1]

# Calculate factor scores for collective efficacy using regression method
collective_efficacy_scores <- factor.scores(collective_efficacy_items, fa_collective_efficacy, method = "regression")$scores[, 1]

# Add factor scores back to the dataset
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    physical_resources = physical_resources_scores,
    environmental_threats = environmental_threats_scores,
    collective_efficacy = collective_efficacy_scores
  )

# Remove intermediate data frames from the global environment
rm(physical_resources_items, environmental_threats_items, collective_efficacy_items, 
   fa_physical_resources, fa_environmental_threats, fa_collective_efficacy, 
   physical_resources_scores, environmental_threats_scores, collective_efficacy_scores)

Child Health & Development Variables

Standardize Down Syndrome & Tourette Syndrome severity variables based on survey year, calculate total severity scores across diagnoses, categorize average severity scores into severity levels, recode child general health, and convert the special healthcare needs indicator into binary format with labeled categories.

# Standardize Down Syndrome and Tourette Syndrome severity variables based on survey year
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    downsyn_desc_recode = case_when(
      survey_year %in% c(2016, 2017) & downsyn_desc == 1 ~ 1,
      survey_year %in% c(2016, 2017) & downsyn_desc == 4 ~ 3,
      survey_year >= 2019 ~ downsyn_desc,
      TRUE ~ NA_real_  # Handle missing or undefined years as NA
    ),
    k2q38c_recode = case_when(
      k2q38c == 1 ~ 1,
      k2q38c == 4 ~ 3
    )
  )

# List of severity variables, including recoded Down & Tourette Syndromes' description
severity_variables <- c("k2q31c", "k2q35c", "k2q38c_recode", "k2q30c", "k2q36c", "k2q60c", 
                        "k2q37c", "downsyn_desc_recode", "cerpals_desc", "k2q42c")

# Calculate total and average severity scores
neurodivergent_data <- neurodivergent_data %>%
  rowwise() %>%
  mutate(
    total_severity_score = ifelse(all(is.na(c_across(all_of(severity_variables)))), NA_real_, 
                                  sum(c_across(all_of(severity_variables)), na.rm = TRUE))
  ) %>%
  ungroup()

# Recode child general health variable
neurodivergent_data <- neurodivergent_data %>%
  mutate(child_gen_health = factor(
    case_when(
      k2q01 == 1 ~ "Excellent",
      k2q01 == 2 ~ "Very Good",
      k2q01 == 3 ~ "Good",
      k2q01 == 4 ~ "Fair",
      k2q01 == 5 ~ "Poor"
    ),
    levels = c("Excellent", "Very Good", "Good", "Fair", "Poor"))
    )

# Convert Special Healthcare Needs (SHCN) indicator to binary and label accordingly
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    sc_cshcn = ifelse(sc_cshcn == 1, 1, 0),
    special_healthcare_needs = factor(
      case_when(
        sc_cshcn == 0 ~ "Non-SHCN",
        sc_cshcn == 1 ~ "SHCN"
      ),
      levels = c("Non-SHCN", "SHCN")
    )
  )

1

k2q31c = ADD/ADHD Severity Description; k2q35c = Autism ASD Severity Description; k2q38c_recode = Tourette Syndrome Severity Description; k2q30c = Learning Disability Severity Description; k2q36c = Developmental Delay Severity Description; k2q60c = Intellectual Disability Severity Description; k2q37c = Speech Disorder Severity Description; downsyn_desc_recode = Down Syndrome Severity Description; cerpals_desc = Cerebral Palsy Severity Description; k2q42c = Epilepsy Severity Description

Psychosocial Risk & Protective Factors

Recode ACEs, psychopathology, and psychosocial factors. Compute the Childhood Flourishing Index and the Family Resilience and Connection Index, based on Bethell et al. (2019), “Family Resilience And Connection Promote Flourishing Among US Children, Even Amid Adversity.” Create categorical versions of these indices for further analysis.

• The Child Flourishing Index (CFI) in was calculated using three items from the National Survey of Children’s Health (NSCH). Each item asked parents how well each of the following statements described their child: (1) “shows interest and curiosity in learning new things,” (2) “works to finish tasks he or she starts,” and (3) “stays calm and in control when faced with a challenge.” Each item was assigned one point for a response of “definitely true.” Children who scored 3 were classified as flourishing.

• The Family Resilience and Connection Index (FRCI) was created using six items from the NSCH. The index included four items on family resilience, which asked how often families: (1) “talk together about what to do,” (2) “work together to solve our problems,” (3) “know we have strengths to draw on,” and (4) “stay hopeful even in difficult times.” Parents could earn one point for each of these items if they responded “all of the time.” Two additional items assessed family connection: (1) the ability of the parent to “share ideas or talk about things that really matter” with their child and (2) how well parents were “handling the day-to-day demands of raising children.” One point was assigned for each response of “very well.” The FRCI score ranged from 0 to 6.

# Recode Adverse Childhood Experiences (ACE) variables & psychopathology outcomes
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Recode 'ace1' to a dichotomous variable
    ace1_recode = ifelse(ace1 %in% c(2, 3, 4), 1, ifelse(ace1 == 1, 0, NA))
  ) %>%
  mutate(
    # Total ACE count excluding missing values
    ace_total = rowSums(select(., ace1_recode, ace3:ace12) == 1, na.rm = TRUE),
    
    # Categorize total ACEs for analysis
    ace_counts = cut(ace_total, breaks = c(-1, 0, 1, 3, Inf), labels = c('0 ACEs', '1 ACE', '2-3 ACEs', '4+ ACEs'), right = TRUE),
    ace_counts = factor(ace_counts, levels = c("0 ACEs", "1 ACE", "2-3 ACEs", "4+ ACEs")),
    
    # Recode psychopathology outcomes
    Anxiety = factor(k2q33b, levels = c(2, 1), labels = c("No", "Yes")),
    Depression = factor(k2q32b, levels = c(2, 1), labels = c("No", "Yes")),
    Behavioral_Problems = factor(k2q34b, levels = c(2, 1), labels = c("No", "Yes"))
  )

# Calculate and recode Childhood Flourishing Index (CFI)
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Calculate the Childhood Flourishing Index (CFI)
    CFI = rowSums(select(., k6q71_r, k7q85_r, k7q84_r) == 1, na.rm = TRUE),
    
    # Recode CFI into a dichotomous variable
    CFI_dich = case_when(
      CFI == 3 ~ "Flourishing",
      CFI %in% 0:2 ~ "Not Flourishing",
      TRUE ~ NA_character_
    ),
    CFI_dich = factor(CFI_dich, levels = c("Not Flourishing", "Flourishing"))
  )

# Calculate and recode Family Resilience and Connection Index (FRCI)
neurodivergent_data <- neurodivergent_data %>%
  mutate(
    # Recode FRCI items: 1 point for "All of the time" (response = 1), handle missing values
    talkabout_score = ifelse(!is.na(talkabout) & talkabout == 1, 1, 0),
    wktosolve_score = ifelse(!is.na(wktosolve) & wktosolve == 1, 1, 0),
    strengths_score = ifelse(!is.na(strengths) & strengths == 1, 1, 0),
    hopeful_score = ifelse(!is.na(hopeful) & hopeful == 1, 1, 0),
    
    # Recode additional FRCI items: 1 point for "Very well" (response = 1), handle missing values
    k8q21_score = ifelse(!is.na(k8q21) & k8q21 == 1, 1, 0),
    k8q30_score = ifelse(!is.na(k8q30) & k8q30 == 1, 1, 0),
    
    # Calculate total FRCI score by summing all items, ignoring missing values
    FRCI = rowSums(across(c(talkabout_score, wktosolve_score, strengths_score, hopeful_score, k8q21_score, k8q30_score)), na.rm = TRUE),
    
    # Create categorical version of FRCI
    FRCI_cat = factor(
      case_when(
        FRCI <= 1 ~ "0-1",
        FRCI <= 3 ~ "2-3",
        FRCI <= 6 ~ "4-6",
        TRUE ~ NA_character_
      ),
      levels = c("0-1", "2-3", "4-6")
    )
  )

1

ace1 is originally coded by NSCH as a frequency measure, but it’s needed as indicative variable of presence like the other ace# variables.

2

ace1 = Hard to Cover Basics Like Food or Housing; ace3 = Child Experienced - Parent or Guardian Divorced; ace4 = Child Experienced - Parent or Guardian Died; ace5 = Child Experienced - Parent or Guardian Time in Jail; ace6 = Child Experienced - Adults Slap, Hit, Kick, Punch Others; ace7 = Child Experienced - Victim of Violence; ace8 = Child Experienced - Lived with Mentally Ill; ace9 = Child Experienced - Lived with Person with Alcohol/Drug Problem; ace10 = Child Experienced - Treated Unfairly Because of Race; ace11 = Child Experienced - Treated Unfairly Because of Health Condition; ace12 = Child Experienced - Treated Unfairly Because of their Sexual Orientation or Gender Identity

3

k2q33b = Anxiety Currently; k2q32b = Depression Currently; k2q34b = Behavior Problems Currently

4

k6q71_r = Show Interest and Curiosity; k7q85_r = Stays Calm and In Control When Challenged; k7q84_r = Works to Finish Tasks Started

5

talkabout = Facing Problems - How Often Talk Together; wktosolve = Facing Problems - How Often Work Together; strengths = Facing Problems - How Often Draw on Strengths; k8q21 = Share Ideas or Talk About Things That Matter; k8q30 = How Well Handling Demands of Raising Children

Frequencies & Descriptive Statistics

---

Sample Trends of Neurodevelopmental Diagnoses

Calculate yearly prevalence rates for neurodevelopmental disorder (NDD) diagnoses without survey weights. Group data by survey year, calculate prevalence for each diagnosis, and transform it into a long format. Create a line plot to visualize yearly trends in the prevalence of each NDD diagnosis, using color to differentiate between the diagnoses.

# Generate list of NDD diagnoses among subject sample
variables <- c("k2q31a", "k2q35a", "k2q38a", "k2q30a", "k2q36a", "k2q60a", "k2q37a", 
               "downsyn", "k2q61a", "k2q42a")
variable_names <- c("ADHD", "Autism/ASD", "Tourette’s Syndrome", "Learning Disability",
                    "Development Delay", "Intellectual Disability", 
                    "Speech/Other Language Disorder", "Down Syndrome", 
                    "Cerebral Palsy", "Epilepsy")

# Calculate yearly prevalence rates without survey weights for specific sample
sample_prevalence <- neurodivergent_data %>%
  group_by(survey_year) %>%
  summarise(across(all_of(variables), 
                   ~ mean(. == 1, na.rm = TRUE) * 100, 
                   .names = "Prevalence_{col}")) %>%
  pivot_longer(
    cols = starts_with("Prevalence_"), 
    names_to = "Diagnosis", 
    values_to = "Prevalence"
  ) %>%
  mutate(
    Diagnosis = gsub("Prevalence_", "", Diagnosis),
    Diagnosis = factor(Diagnosis, levels = variables, labels = variable_names)
  )

# Plot the prevalence rates by year for each NDD diagnosis without survey weights
ggplot(sample_prevalence, aes(x = survey_year, y = Prevalence, group = Diagnosis, color = Diagnosis)) +
  geom_line(size = 1.1) +
  geom_point(size = 2) +
  labs(
    title = "Yearly Prevalence of Neurodevelopmental Disorders (Sample Estimates)",
    x = "Survey Year",
    y = "Prevalence (%)",
    color = "Diagnosis"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

1

The following diagnoses’ variables correspond to whether a child has EVER received such a diagnosis: k2q31a = ADD/ADHD; k2q35a = Autism ASD; k2q38a = Tourette Syndrome; k2q30a = Learning Disability; k2q36a = Developmental Delay; k2q60a = Intellectual Disability; k2q37a = Speech Disorder; downsyn = Down Syndrome; cerpals = Cerebral Palsy; k2q42a = Epilepsy

# Remove intermediate data frames from the global environment
rm(design, sample_prevalence, yearly_prevalence, yearly_prevalence_df)

Sample Demographic Distributions

Calculate and visualize the distribution of child age, educational attainment, sex, socioeconomic status, and race among neurodivergent participants. Create a raincloud plot for age distribution alongside a summary statistics table, and generate bar plots for key categorical variables, ensuring all visuals are clear and informative.

# Calculate descriptives of age for neurodivergent participants
neurodiv_age_summary <- neurodivergent_data %>%
  summarise(
    Mean = mean(sc_age_years, na.rm = TRUE),
    n = n(),
    SE = sd(sc_age_years, na.rm = TRUE) / sqrt(n()),
    SD = sd(sc_age_years, na.rm = TRUE),
    Minimum = min(sc_age_years, na.rm = TRUE),
    Q1 = quantile(sc_age_years, 0.25, na.rm = TRUE),
    Median = median(sc_age_years, na.rm = TRUE),
    Q3 = quantile(sc_age_years, 0.75, na.rm = TRUE),
    Maximum = max(sc_age_years, na.rm = TRUE),
    IQR = quantile(sc_age_years, 0.75, na.rm = TRUE) - quantile(sc_age_years, 0.25, na.rm = TRUE),
    .groups = 'drop'
  )

# Create raincloud plot displaying child age descriptives
raincloud_plot <- neurodivergent_data %>%
  ggplot(aes(x = 1, y = sc_age_years)) +
  PupillometryR::geom_flat_violin(
    adjust = 1.5, trim = FALSE, alpha = 0.5, color = NA, fill = "deeppink", 
    position = position_nudge(x = 0.2)
  ) +
  geom_boxplot(
    width = 0.12, outlier.color = "black", colour = "black", fill = NA, alpha = 0.5, 
    position = position_nudge(x = 0.2)
  ) +
  stat_dots(
    dotsize = 0.1, side = "left", justification = 1, binwidth = 1, alpha = 0.5, 
    colour = "deeppink", fill = "deeppink"
  ) +
  geom_point(data = neurodiv_age_summary, aes(x = 1, y = Mean), shape = 18, 
             color = "deeppink", size = 3, position = position_nudge(x = 0.2)) +
  geom_errorbar(data = neurodiv_age_summary, aes(x = 1, y = Mean, ymin = Mean - SE, ymax = Mean + SE), 
                width = 0.1, color = "deeppink", position = position_nudge(x = 0.2)) +
  labs(title = "Age Distribution of Neurodivergent Participants", y = "Age (Years)", x = "") +
  scale_y_continuous(breaks = seq(5, 18, by = 1), limits = c(5, 18)) +
  coord_flip() +
  theme_minimal() +
  theme(axis.text.y = element_blank(), legend.position = "none")

# Create summary statistics table as a data frame and transpose it
summary_df <- neurodiv_age_summary %>%
  select(Mean, SD, SE, Median, Minimum, Maximum, IQR, Q1, Q3) %>%
  t() %>%
  as.data.frame() %>%
  mutate_all(~ round(., 2))

# Update row names to meaningful labels
colnames(summary_df) <- "Value"
summary_df$Statistic <- rownames(summary_df)
summary_df <- summary_df %>%
  select(Statistic, Value)

# Convert the transposed summary statistics data frame to a table using gridExtra
summary_table <- gridExtra::tableGrob(summary_df, rows = NULL)

# Combine the plot and summary table using gridExtra
gridExtra::grid.arrange(raincloud_plot, summary_table, nrow = 1)

# Bar plot for educational attainment counts categories
ggplot(neurodivergent_data %>% filter(!is.na(caregiver_education)), aes(x = caregiver_education)) +
  geom_bar(aes(fill = caregiver_education), show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -0.5) +
  labs(
    title = "Distribution of Adults' Highest Educational Attainment",
    x = "Educational Attainment",
    y = "Count"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Bar plot for sex counts distribution
ggplot(neurodivergent_data, aes(x = sc_sex_cat)) +
  geom_bar(aes(fill = sc_sex_cat), show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  scale_fill_manual(values = c("Male" = "lightblue", "Female" = "pink")) +
  labs(
    title = "Sex Distribution of Neurodivergent Participants",
    x = "Sex",
    y = "Count"
  ) +
  theme_minimal()

# Bar plot for SES counts categories
ggplot(neurodivergent_data, aes(x = fpl_category)) +
  geom_bar(aes(fill = fpl_category), show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  labs(
    title = "Distribution of Socioeconomic Status (FPL)",
    x = "Federal Poverty Level Groups",
    y = "Count"
  ) +
  theme_minimal()

# Bar plot for race counts categories
ggplot(neurodivergent_data, aes(x = sc_race_cat)) +
  geom_bar(aes(fill = sc_race_cat), show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), vjust = -0.3, size = 3.5) +
  labs(
    title = "Race Distribution of Neurodivergent Participants",
    x = "Race",
    y = "Count"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Remove intermediate data frames from the global environment
rm(neurodiv_age_summary, raincloud_plot, summary_df, summary_table)

Sample Trends of Child Flourishing

Create visualizations for the Childhood Flourishing Index (CFI) and dichotomous CFI variable (Flourishing vs Not Flourishing). Generate bar plots to show the distribution of these variables and line plots to display yearly prevalence rates across survey years.

# Bar plot for Childhood Flourishing Index (CFI) Score distribution
ggplot(neurodivergent_data, aes(x = factor(CFI), fill = factor(CFI))) +
  geom_bar(show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  labs(
    x = "Childhood Flourishing Index (CFI) Score",
    y = "Number of Children",
    title = "Distribution of Children by CFI Score"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0))

# Calculate prevalence of CFI scores grouped by year
flourishing_prevalence <- neurodivergent_data %>%
  group_by(survey_year, CFI) %>%
  summarise(
    Count = n(),
    .groups = "drop"
  ) %>%
  group_by(survey_year) %>%
  mutate(
    Total = sum(Count),
    Prevalence = (Count / Total) * 100
  )

# Line plot for yearly prevalence of CFI scores
ggplot(flourishing_prevalence, aes(x = survey_year, y = Prevalence, color = as.factor(CFI), group = CFI)) +
  geom_line(size = 1.2) +
  geom_point(size = 3) +
  labs(
    title = "Yearly Prevalence of Child Flourishing Index Scores (Sample Estimate)",
    x = "Survey Year",
    y = "Prevalence (%)",
    color = "Child Flourishing Index Score"
  ) +
  theme_minimal()

# Bar plot for CFI dichotomous variable distribution
ggplot(neurodivergent_data, aes(x = CFI_dich, fill = CFI_dich)) +
  geom_bar(show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  labs(
    title = "Distribution of CFI Dichotomous Variable",
    x = "Child Flourishing Index (Dichotomous)",
    y = "Count"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0))

# Calculate prevalence of flourishing status (CFI_dich) grouped by year
flourishing_prevalence_status <- neurodivergent_data %>%
  group_by(survey_year, CFI_dich) %>%
  summarise(
    Count = n(),
    .groups = "drop"
  ) %>%
  group_by(survey_year) %>%
  mutate(
    Total = sum(Count),
    Prevalence = (Count / Total) * 100
  )

# Line plot for yearly prevalence of flourishing status
ggplot(flourishing_prevalence_status, aes(x = survey_year, y = Prevalence, color = CFI_dich, group = CFI_dich)) +
  geom_line(size = 1.2) +
  geom_point(size = 3) +
  labs(
    title = "Yearly Prevalence of Child Flourishing Status (Sample Estimate)",
    x = "Survey Year",
    y = "Prevalence (%)",
    color = "Child Flourishing Status"
  ) +
  theme_minimal()

# Remove intermediate data frames from the global environment
rm(flourishing_prevalence, flourishing_prevalence_status)

Sample Trends of Family Resilience & Connection

Create visualizations for the Family Resilience and Connection Index (FRCI) scores and categories. Generate bar plots to show the distribution of FRCI scores and categories, and line plots to display yearly prevalence rates across survey years for both continuous FRCI scores and categorical FRCI status.

# Bar plot for Family Resilience and Connection Index (FRCI) Score distribution
ggplot(neurodivergent_data, aes(x = factor(FRCI), fill = factor(FRCI))) +
  geom_bar(show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  labs(
    x = "Family Resilience and Connection Index (FRCI) Score",
    y = "Number of Children",
    title = "Distribution of Children by FRCI Score"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0))

# Calculate prevalence of FRCI scores grouped by year
frci_prevalence_scores <- neurodivergent_data %>%
  group_by(survey_year, FRCI) %>%
  summarise(
    Count = n(),
    .groups = "drop"
  ) %>%
  group_by(survey_year) %>%
  mutate(
    Total = sum(Count),
    Prevalence = (Count / Total) * 100
  )

# Line plot for yearly prevalence of FRCI scores
ggplot(frci_prevalence_scores, aes(x = survey_year, y = Prevalence, color = as.factor(FRCI), group = FRCI)) +
  geom_line(size = 1.2) +
  geom_point(size = 3) +
  labs(
    title = "Yearly Prevalence of Family Resilience & Connection Index (FRCI) Scores (Sample Estimate)",
    x = "Survey Year",
    y = "Prevalence (%)",
    color = "FRCI Score"
  ) +
  theme_minimal()

# Bar plot for FRCI categorical variable distribution
ggplot(neurodivergent_data, aes(x = FRCI_cat, fill = FRCI_cat)) +
  geom_bar(show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  labs(
    title = "Distribution of FRCI Categorical Variable",
    x = "Family Resilience and Connection Index (Categorical)",
    y = "Count"
  ) +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0))

# Calculate prevalence of FRCI categories grouped by year
frci_prevalence_categories <- neurodivergent_data %>%
  group_by(survey_year, FRCI_cat) %>%
  summarise(
    Count = n(),
    .groups = "drop"
  ) %>%
  group_by(survey_year) %>%
  mutate(
    Total = sum(Count),
    Prevalence = (Count / Total) * 100
  )

# Line plot for yearly prevalence of FRCI categories
ggplot(frci_prevalence_categories, aes(x = survey_year, y = Prevalence, color = FRCI_cat, group = FRCI_cat)) +
  geom_line(size = 1.2) +
  geom_point(size = 3) +
  labs(
    title = "Yearly Prevalence of Family Resilience & Connection Index Categories (Sample Estimate)",
    x = "Survey Year",
    y = "Prevalence (%)",
    color = "FRCI Category"
  ) +
  theme_minimal()

# Remove intermediate data frames from the global environment
rm(frci_prevalence_scores, frci_prevalence_categories)

Sample Trends of Psychopathology

Create a stacked bar plot to show the distribution of psychopathology outcomes across participants. Calculate the yearly prevalence of anxiety, depression, and behavioral problems, and visualize these trends using a line plot.

# Create a long-format data frame for plotting psychopathology outcomes
psychopathology_data <- neurodivergent_data %>%
  select(Anxiety, Depression, Behavioral_Problems) %>%
  pivot_longer(cols = everything(), names_to = "Psychopathology", values_to = "Response") %>%
  filter(!is.na(Response)) # Exclude missing responses

# Create a stacked bar plot for psychopathology outcome distribution
ggplot(psychopathology_data, aes(x = Psychopathology, fill = Response)) +
  geom_bar(position = "stack") +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  scale_fill_manual(values = c("Yes" = "blue", "No" = "red")) +
  labs(
    title = "Distribution of Psychopathology Outcomes",
    x = "Psychopathology Outcome",
    y = "Count",
    fill = "Response"
  ) +
  theme_minimal()

# List of psychopathology outcomes for analysis
psychopathology_variables <- c("Anxiety", "Depression", "Behavioral_Problems")

# Calculate prevalence of each psychopathology outcome by survey year
psychopathology_prevalence <- lapply(psychopathology_variables, function(var) {
  neurodivergent_data %>%
    group_by(survey_year, !!sym(var)) %>%
    summarise(
      Count = n(), # Count of each level
      .groups = "drop"
    ) %>%
    group_by(survey_year) %>%
    mutate(
      Total = sum(Count), # Total count per year
      Prevalence = (Count / Total) * 100, # Calculate prevalence as a percentage
      Outcome = var # Add the outcome name
    ) %>%
    filter(!!sym(var) == "Yes") # Only consider "Yes" (those with the outcome)
}) %>%
  bind_rows() # Combine all results into one data frame

# Visualize the prevalence of psychopathology outcomes by year using a line plot
ggplot(psychopathology_prevalence, aes(x = survey_year, y = Prevalence, color = Outcome, group = Outcome)) +
  geom_line(size = 1.2) + # Add lines
  geom_point(size = 3) + # Add points
  labs(
    title = "Yearly Prevalence of Psychopathology Outcomes (Sample Estimate)",
    x = "Survey Year",
    y = "Prevalence (%)",
    color = "Outcome"
  ) +
  theme_minimal() +
  theme(
    axis.text.x = element_text(angle = 0)
  )

# Remove intermediate data frames from the global environment
rm(psychopathology_data, psychopathology_prevalence)

Distribution & Cross-Tabulations of ACEs

Create a bar plot to visualize the distribution of children by ACE counts. Conduct cross-tabulations between ACE counts and psychopathology/flourishing variables, as well as between psychopathology outcomes and the flourishing variable, displaying the results in table format.

# Bar plot of ACE counts categories
ggplot(neurodivergent_data, aes(x = factor(ace_counts), fill = factor(ace_counts))) +
  geom_bar(show.legend = FALSE) +
  geom_text(stat = 'count', aes(label = ..count..), position = position_stack(vjust = 0.5), color = "white") +
  labs(x = "ACE Counts", y = "Number of Children", fill = "ACE Counts",
       title = "Distribution of Children by ACE Counts") +
  theme_minimal() +
  theme(axis.text.x = element_text(angle = 0))

1

Cross-tabulations (or contingency tables) summarize the relationship between two or more categorical variables. In this case, I’m exploring how the mental health conditions (Anxiety, Depression, Behavioral Issues) and Child Flourishing Index (CFI) are distributed across different levels of ACEs.

# Cross-tabulation between psychopathology/flourishing variables and ACEs count variable
crosstable(neurodivergent_data, c(Anxiety, Depression, Behavioral_Problems, CFI_dich), by = ace_counts, showNA = "no", percent_digits = 1, total = "both") %>%
  as_flextable()

# Cross-tabulation between psychopathology outcomes and flourishing variable
crosstable(neurodivergent_data, c(Anxiety, Depression, Behavioral_Problems), by = CFI_dich, showNA = "no", percent_digits = 1, total = "both") %>%
  as_flextable()

label	variable	ace_counts				Total
		0 ACEs	1 ACE	2-3 ACEs	4+ ACEs
Anxiety	No	383 (27.0%)	375 (26.5%)	407 (28.7%)	252 (17.8%)	1417 (7.7%)
	Yes	3389 (20.0%)	4366 (25.8%)	5180 (30.6%)	3992 (23.6%)	16927 (92.3%)
	Total	3772 (20.6%)	4741 (25.8%)	5587 (30.5%)	4244 (23.1%)	18344 (100.0%)
Depression	No	310 (20.3%)	349 (22.8%)	486 (31.8%)	383 (25.1%)	1528 (16.6%)
	Yes	1004 (13.1%)	1475 (19.3%)	2553 (33.4%)	2619 (34.2%)	7651 (83.4%)
	Total	1314 (14.3%)	1824 (19.9%)	3039 (33.1%)	3002 (32.7%)	9179 (100.0%)
Behavioral_Problems	No	592 (21.7%)	692 (25.4%)	875 (32.1%)	569 (20.9%)	2728 (14.5%)
	Yes	3052 (18.9%)	4234 (26.3%)	4939 (30.6%)	3890 (24.1%)	16115 (85.5%)
	Total	3644 (19.3%)	4926 (26.1%)	5814 (30.9%)	4459 (23.7%)	18843 (100.0%)
CFI_dich	Not Flourishing	13914 (27.4%)	14990 (29.6%)	13868 (27.3%)	7936 (15.7%)	50708 (93.0%)
	Flourishing	1530 (40.1%)	1239 (32.4%)	808 (21.2%)	243 (6.4%)	3820 (7.0%)
	Total	15444 (28.3%)	16229 (29.8%)	14676 (26.9%)	8179 (15.0%)	54528 (100.0%)

label	variable	CFI_dich		Total
		Not Flourishing	Flourishing
Anxiety	No	1332 (94.0%)	85 (6.0%)	1417 (7.7%)
	Yes	16579 (97.9%)	348 (2.1%)	16927 (92.3%)
	Total	17911 (97.6%)	433 (2.4%)	18344 (100.0%)
Depression	No	1477 (96.7%)	51 (3.3%)	1528 (16.6%)
	Yes	7530 (98.4%)	121 (1.6%)	7651 (83.4%)
	Total	9007 (98.1%)	172 (1.9%)	9179 (100.0%)
Behavioral_Problems	No	2633 (96.5%)	95 (3.5%)	2728 (14.5%)
	Yes	15943 (98.9%)	172 (1.1%)	16115 (85.5%)
	Total	18576 (98.6%)	267 (1.4%)	18843 (100.0%)

Inferential Statistics

---

Interaction Models

Fit three logistic regression models with survey weights to examine the moderating effect of Child Flourishing Index (CFI) on the association between Adverse Childhood Experiences (ACEs) and psychopathology outcomes (anxiety, depression, behavioral problems). Extract residuals, check heteroscedasticity using Breusch-Pagan tests, calculate Cook’s distance for influential observations, and obtain robust standard errors. Create summary tables with robust standard errors and False Discovery Rate (FDR) adjusted p-values for each model.

Anxiety Model

# Create survey design object with Final Weight for Child (FWC) as weights
design <- svydesign(ids = ~1, weights = ~fwc, data = neurodivergent_data) # will be used across all models

# Logistic regression model with CFI_dich as a moderator for Anxiety, using robust standard errors
model_anx_mod <- svyglm(Anxiety ~ ace_counts * CFI_dich + sc_race_cat + sc_sex_cat + 
                        fpl_category + caregiver_education + sc_age_years + survey_year + 
                        caregiver_mental_health + child_gen_health + 
                        caregiver_composition + physical_resources + 
                        environmental_threats + collective_efficacy +
                        special_healthcare_needs + total_severity_score + 
                        HAI_category + barriers_cat + insurance_type + IAI_category, 
                        design = design, family = quasibinomial())

# Extract residuals from the logistic regression model
residuals_anx <- residuals(model_anx_mod, type = "deviance")

# Plot residuals to visually inspect them
plot(residuals_anx, main = "Residuals Plot for Anxiety Model", ylab = "Residuals", xlab = "Index", pch = 20, col = "blue")
abline(h = 0, col = "red", lty = 2)

# Perform the Breusch-Pagan test for heteroscedasticity
bptest(model_anx_mod)

# Calculate Cook's distance for influential observations (leverage) in each model
cooks_anx <- cooks.distance(model_anx_mod)

# Plot for Anxiety model
plot(cooks_anx, type = "p", main = "Cook's Distance for Anxiety Model",
     ylab = "Cook's Distance", xlab = "Observation Index")
abline(h = 4 / length(cooks_anx), col = "red", lty = 2) # Add a reference line for influential points

# Obtain robust standard errors for Anxiety model
robust_anx <- coeftest(model_anx_mod, vcov = vcovHC(model_anx_mod, type = "HC0"))

# Create a summary table for the Anxiety model with robust standard errors
tbl_anx_mod <- tbl_regression(model_anx_mod, exponentiate = TRUE, label = list(
      ace_counts ~ "ACEs",
      CFI_dich ~ "Child Flourishing Index",
      fpl_category ~ "Federal Poverty Level",
      caregiver_education ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex",
      survey_year ~ "Survey Year",
      caregiver_mental_health ~ "Caregiver Mental Health",
      child_gen_health ~ "Child General Health",
      caregiver_composition ~ "Caregiver Composition",
      physical_resources ~ "Neighborhood Physical Resources",
      environmental_threats ~ "Neighborhood Environmental Threats",
      collective_efficacy ~ "Neighborhood Collective Efficacy",
      special_healthcare_needs ~ "Special Healthcare Needs",
      total_severity_score ~ "NDD Severity (Total)",
      HAI_category ~ "Healthcare Access Index",
      barriers_cat ~ "Access Barriers Index",
      insurance_type ~ "Insurance Type",
      IAI_category ~ "Insurance Adequacy Index"
    )) %>%
  add_q() %>% # adjusts p-values for multiple testing (FDR)
  bold_p(t = 0.05) %>%
  bold_p(t = 0.05, q = TRUE) %>% # now bold q-values under the threshold of 0.05
  italicize_levels() %>%
  gtsummary::as_gt() %>%
  gt::tab_source_note(gt::md("*Note*. This logistic regression model examines the moderating effect of the Child Flourishing Index (CFI) on the association between Adverse Childhood Experiences (ACEs) and anxiety outcomes in neurodivergent children. Survey weights from the National Survey of Children's Health (NSCH) were applied to ensure the representativeness of the target population. Robust standard errors (HC0) were used to account for potential heteroscedasticity not fully addressed by the survey weights. HC0 provides the original White standard errors, which are suitable for large samples and do not include small-sample bias corrections. Covariates include child demographics, family characteristics, neighborhood conditions and survey year, with FDR-adjusted p-values to account for multiple comparisons."))

# Print table to review results
tbl_anx_mod

    studentized Breusch-Pagan test

data:  model_anx_mod
BP = 851.64, df = 52, p-value < 2.2e-16

Characteristic	OR1	95% CI1	p-value	q-value2
ACEs
0 ACEs	—	—
1 ACE	1.14	0.76, 1.71	0.5	0.8
2-3 ACEs	1.40	0.89, 2.18	0.14	0.3
4+ ACEs	1.55	0.95, 2.54	0.079	0.2
Child Flourishing Index
Not Flourishing	—	—
Flourishing	1.18	0.48, 2.89	0.7	0.9
Child Race
Non-Hispanic White	—	—
Non-Hispanic Black or African American	0.78	0.45, 1.36	0.4	0.6
Non-Hispanic American Indian or Alaska Native	0.57	0.23, 1.44	0.2	0.5
Non-Hispanic Asian, Native Hawaiian, or Pacific Islander	0.33	0.17, 0.65	0.001	0.012
Non-Hispanic Other or Mixed Race	0.94	0.61, 1.46	0.8	0.9
Hispanic/Latino of Any Race	0.61	0.43, 0.86	0.005	0.023
Child Sex
Male	—	—
Female	0.97	0.72, 1.31	0.9	0.9
Federal Poverty Level
<100%	—	—
100%-199%	0.58	0.34, 1.00	0.049	0.15
200%-299%	0.76	0.42, 1.36	0.4	0.6
300%-399%	0.51	0.28, 0.93	0.028	0.10
≥400%	0.75	0.42, 1.34	0.3	0.6
Caregiver Education
Less Than High School	—	—
High School/Vocational	3.11	1.42, 6.80	0.004	0.021
Some College/Associate Degree	2.99	1.44, 6.21	0.003	0.021
College Degree/Higher	4.25	2.00, 9.01	<0.001	0.003
Child Age	0.99	0.95, 1.03	0.7	0.9
Survey Year
2016	—	—
2017	1.14	0.61, 2.14	0.7	0.9
2018	0.92	0.56, 1.52	0.8	0.9
2019	0.65	0.40, 1.07	0.087	0.2
2020	1.21	0.68, 2.15	0.5	0.8
2021	1.90	1.22, 2.96	0.004	0.021
2022	1.44	0.89, 2.33	0.14	0.3
2023	1.24	0.76, 2.02	0.4	0.6
Caregiver Mental Health
Excellent	—	—
Very Good	1.25	0.90, 1.73	0.2	0.4
Good	1.78	1.21, 2.60	0.003	0.021
Fair	2.45	1.43, 4.20	0.001	0.012
Poor	5.86	1.74, 19.7	0.004	0.021
Child General Health
Excellent	—	—
Very Good	1.03	0.77, 1.38	0.9	0.9
Good	1.58	1.05, 2.39	0.029	0.10
Fair	2.02	0.86, 4.76	0.11	0.3
Poor	1.55	0.41, 5.79	0.5	0.8
Caregiver Composition
Dual-caregiver household	—	—
Single-caregiver household	1.00	0.70, 1.43	>0.9	>0.9
Neighborhood Physical Resources	0.81	0.69, 0.96	0.012	0.049
Neighborhood Environmental Threats	1.01	0.85, 1.20	>0.9	>0.9
Neighborhood Collective Efficacy	0.86	0.73, 1.02	0.086	0.2
Special Healthcare Needs
Non-SHCN	—	—
SHCN	1.78	1.25, 2.53	0.001	0.012
NDD Severity (Total)	1.15	1.08, 1.21	<0.001	<0.001
Healthcare Access Index
High Access	—	—
Medium Access	0.36	0.25, 0.51	<0.001	<0.001
Low Access	0.92	0.53, 1.59	0.8	0.9
Access Barriers Index
No Barriers	—	—
Few Barriers	0.87	0.56, 1.37	0.6	0.8
Moderate Barriers	0.89	0.38, 2.08	0.8	0.9
Many Barriers	0.75	0.19, 2.92	0.7	0.9
Insurance Type
Private Only	—	—
Private + Public	1.06	0.61, 1.84	0.8	0.9
Public Only	0.70	0.47, 1.05	0.081	0.2
Uninsured	0.09	0.01, 0.94	0.044	0.14
Insurance Adequacy Index
High Adequacy	—	—
Moderate Adequacy	1.22	0.85, 1.76	0.3	0.5
Low Adequacy	0.91	0.57, 1.46	0.7	0.9
ACEs * Child Flourishing Index
1 ACE * Flourishing	1.93	0.52, 7.23	0.3	0.6
2-3 ACEs * Flourishing	0.80	0.16, 3.92	0.8	0.9
4+ ACEs * Flourishing	1.74	0.26, 11.6	0.6	0.8
Note. This logistic regression model examines the moderating effect of the Child Flourishing Index (CFI) on the association between Adverse Childhood Experiences (ACEs) and anxiety outcomes in neurodivergent children. Survey weights from the National Survey of Children’s Health (NSCH) were applied to ensure the representativeness of the target population. Robust standard errors (HC0) were used to account for potential heteroscedasticity not fully addressed by the survey weights. HC0 provides the original White standard errors, which are suitable for large samples and do not include small-sample bias corrections. Covariates include child demographics, family characteristics, neighborhood conditions and survey year, with FDR-adjusted p-values to account for multiple comparisons.
1 OR = Odds Ratio, CI = Confidence Interval
2 False discovery rate correction for multiple testing

Depression Model

# Logistic regression model with CFI_dich as a moderator for Depression, using robust standard errors
model_dep_mod <- svyglm(Depression ~ ace_counts * CFI_dich + sc_race_cat + sc_sex_cat + 
                        fpl_category + caregiver_education + sc_age_years + survey_year + 
                        caregiver_mental_health + child_gen_health + 
                        caregiver_composition + physical_resources + 
                        environmental_threats + collective_efficacy +
                        special_healthcare_needs + total_severity_score + 
                        HAI_category + barriers_cat + insurance_type + IAI_category,  
                        design = design, family = quasibinomial())

# Extract residuals from the logistic regression model
residuals_dep <- residuals(model_dep_mod, type = "deviance")

# Plot residuals to visually inspect them
plot(residuals_dep, main = "Residuals Plot for Depression Model", ylab = "Residuals", xlab = "Index", pch = 20, col = "blue")
abline(h = 0, col = "red", lty = 2)

# Perform the Breusch-Pagan test for heteroscedasticity
bptest(model_dep_mod)

# Calculate Cook's distance for influential observations (leverage) in each model
cooks_dep <- cooks.distance(model_dep_mod)

# Plot for Depression model
plot(cooks_dep, type = "p", main = "Cook's Distance for Depression Model",
     ylab = "Cook's Distance", xlab = "Observation Index")
abline(h = 4 / length(cooks_dep), col = "red", lty = 2) # Add a reference line for influential points

# Obtain robust standard errors for Depression model
robust_dep <- coeftest(model_dep_mod, vcov = vcovHC(model_dep_mod, type = "HC0"))

# Create a summary table for the Depression model with robust standard errors
tbl_dep_mod <- tbl_regression(model_dep_mod, exponentiate = TRUE, label = list(
      ace_counts ~ "ACEs",
      CFI_dich ~ "Child Flourishing Index",
      fpl_category ~ "Federal Poverty Level",
      caregiver_education ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex",
      survey_year ~ "Survey Year",
      caregiver_mental_health ~ "Caregiver Mental Health",
      child_gen_health ~ "Child General Health",
      caregiver_composition ~ "Caregiver Composition",
      physical_resources ~ "Neighborhood Physical Resources",
      environmental_threats ~ "Neighborhood Environmental Threats",
      collective_efficacy ~ "Neighborhood Collective Efficacy",
      special_healthcare_needs ~ "Special Healthcare Needs",
      total_severity_score ~ "NDD Severity (Total)",
      HAI_category ~ "Healthcare Access Index",
      barriers_cat ~ "Access Barriers Index",
      insurance_type ~ "Insurance Type",
      IAI_category ~ "Insurance Adequacy Index"
    )) %>%
  add_q() %>% # adjusts p-values for multiple testing (FDR)
  bold_p(t = 0.05) %>%
  bold_p(t = 0.05, q = TRUE) %>% # now bold q-values under the threshold of 0.05
  italicize_levels() %>%
  gtsummary::as_gt() %>%
  gt::tab_source_note(gt::md("*Note*. This logistic regression model examines the moderating effect of the Child Flourishing Index (CFI) on the association between Adverse Childhood Experiences (ACEs) and depression outcomes in neurodivergent children. Survey weights from the National Survey of Children's Health (NSCH) were applied to ensure the representativeness of the target population. Robust standard errors (HC0) were used to account for potential heteroscedasticity not fully addressed by the survey weights. HC0 provides the original White standard errors, which are suitable for large samples and do not include small-sample bias corrections. Covariates include child demographics, family characteristics, neighborhood conditions and survey year, with FDR-adjusted p-values to account for multiple comparisons."))

# Print table to review results
tbl_dep_mod

    studentized Breusch-Pagan test

data:  model_dep_mod
BP = 468.18, df = 52, p-value < 2.2e-16

Characteristic	OR1	95% CI1	p-value	q-value2
ACEs
0 ACEs	—	—
1 ACE	1.51	1.01, 2.25	0.044	0.14
2-3 ACEs	1.81	1.19, 2.74	0.005	0.039
4+ ACEs	2.15	1.39, 3.30	<0.001	0.009
Child Flourishing Index
Not Flourishing	—	—
Flourishing	1.24	0.37, 4.18	0.7	0.8
Child Race
Non-Hispanic White	—	—
Non-Hispanic Black or African American	0.76	0.46, 1.25	0.3	0.6
Non-Hispanic American Indian or Alaska Native	0.77	0.26, 2.28	0.6	0.8
Non-Hispanic Asian, Native Hawaiian, or Pacific Islander	1.29	0.63, 2.65	0.5	0.8
Non-Hispanic Other or Mixed Race	1.46	0.95, 2.24	0.081	0.2
Hispanic/Latino of Any Race	0.58	0.38, 0.88	0.010	0.058
Child Sex
Male	—	—
Female	1.43	1.07, 1.93	0.017	0.080
Federal Poverty Level
<100%	—	—
100%-199%	1.18	0.72, 1.93	0.5	0.8
200%-299%	0.80	0.44, 1.47	0.5	0.8
300%-399%	0.86	0.42, 1.75	0.7	0.8
≥400%	0.99	0.53, 1.86	>0.9	>0.9
Caregiver Education
Less Than High School	—	—
High School/Vocational	1.14	0.53, 2.41	0.7	0.8
Some College/Associate Degree	0.87	0.43, 1.79	0.7	0.8
College Degree/Higher	0.91	0.41, 2.04	0.8	0.9
Child Age	1.00	0.95, 1.05	>0.9	>0.9
Survey Year
2016	—	—
2017	1.44	0.70, 2.97	0.3	0.6
2018	1.90	1.03, 3.51	0.039	0.14
2019	1.13	0.62, 2.05	0.7	0.8
2020	1.93	1.02, 3.66	0.044	0.14
2021	1.48	0.79, 2.76	0.2	0.5
2022	1.43	0.81, 2.53	0.2	0.5
2023	1.34	0.76, 2.36	0.3	0.6
Caregiver Mental Health
Excellent	—	—
Very Good	1.05	0.74, 1.48	0.8	0.8
Good	1.94	1.30, 2.90	0.001	0.012
Fair	2.08	1.08, 3.99	0.028	0.11
Poor	3.86	1.22, 12.2	0.021	0.093
Child General Health
Excellent	—	—
Very Good	0.92	0.68, 1.25	0.6	0.8
Good	1.71	1.14, 2.57	0.009	0.058
Fair	1.57	0.78, 3.15	0.2	0.5
Poor	8.45	1.98, 36.0	0.004	0.034
Caregiver Composition
Dual-caregiver household	—	—
Single-caregiver household	0.59	0.43, 0.81	0.001	0.012
Neighborhood Physical Resources	0.93	0.78, 1.12	0.5	0.8
Neighborhood Environmental Threats	0.94	0.76, 1.15	0.5	0.8
Neighborhood Collective Efficacy	0.97	0.82, 1.15	0.7	0.8
Special Healthcare Needs
Non-SHCN	—	—
SHCN	2.16	1.48, 3.14	<0.001	0.002
NDD Severity (Total)	1.08	1.02, 1.16	0.016	0.080
Healthcare Access Index
High Access	—	—
Medium Access	0.32	0.22, 0.47	<0.001	<0.001
Low Access	0.73	0.47, 1.12	0.2	0.4
Access Barriers Index
No Barriers	—	—
Few Barriers	0.84	0.60, 1.18	0.3	0.6
Moderate Barriers	0.78	0.34, 1.76	0.5	0.8
Many Barriers	0.81	0.21, 3.07	0.8	0.8
Insurance Type
Private Only	—	—
Private + Public	1.25	0.71, 2.20	0.4	0.8
Public Only	1.07	0.71, 1.61	0.7	0.8
Uninsured	0.22	0.03, 1.57	0.13	0.4
Insurance Adequacy Index
High Adequacy	—	—
Moderate Adequacy	1.17	0.87, 1.58	0.3	0.6
Low Adequacy	1.25	0.76, 2.05	0.4	0.7
ACEs * Child Flourishing Index
1 ACE * Flourishing	0.69	0.11, 4.38	0.7	0.8
2-3 ACEs * Flourishing	4.35	0.88, 21.4	0.070	0.2
4+ ACEs * Flourishing	1.50	0.30, 7.44	0.6	0.8
Note. This logistic regression model examines the moderating effect of the Child Flourishing Index (CFI) on the association between Adverse Childhood Experiences (ACEs) and depression outcomes in neurodivergent children. Survey weights from the National Survey of Children’s Health (NSCH) were applied to ensure the representativeness of the target population. Robust standard errors (HC0) were used to account for potential heteroscedasticity not fully addressed by the survey weights. HC0 provides the original White standard errors, which are suitable for large samples and do not include small-sample bias corrections. Covariates include child demographics, family characteristics, neighborhood conditions and survey year, with FDR-adjusted p-values to account for multiple comparisons.
1 OR = Odds Ratio, CI = Confidence Interval
2 False discovery rate correction for multiple testing

Behavioral Problems Model

# Logistic regression model with CFI_dich as a moderator for Behavioral Problems, using robust standard errors
model_beh_mod <- svyglm(Behavioral_Problems ~ ace_counts * CFI_dich + sc_race_cat + 
                        sc_sex_cat + fpl_category + caregiver_education + sc_age_years + 
                        survey_year + caregiver_mental_health + child_gen_health + 
                        caregiver_composition + physical_resources + environmental_threats +
                        collective_efficacy + special_healthcare_needs + 
                        total_severity_score + HAI_category + barriers_cat + 
                        insurance_type + IAI_category,
                        design = design, family = quasibinomial())

# Extract residuals from the logistic regression model
residuals_beh <- residuals(model_beh_mod, type = "deviance")

# Plot residuals to visually inspect them
plot(residuals_beh, main = "Residuals Plot for Behavioral Problems Model", ylab = "Residuals", xlab = "Index", pch = 20, col = "blue")
abline(h = 0, col = "red", lty = 2)

# Perform the Breusch-Pagan test for heteroscedasticity
bptest(model_beh_mod)

# Calculate Cook's distance for influential observations (leverage) in each model
cooks_beh <- cooks.distance(model_beh_mod)

# Plot for Behavioral Problems model
plot(cooks_beh, type = "p", main = "Cook's Distance for Behavioral Problems Model",
     ylab = "Cook's Distance", xlab = "Observation Index")
abline(h = 4 / length(cooks_beh), col = "red", lty = 2) # Add a reference line for influential points

# Obtain robust standard errors for Behavioral Problems model
robust_beh <- coeftest(model_beh_mod, vcov = vcovHC(model_beh_mod, type = "HC0"))

# Create a summary table for the Behavioral Problems model with robust standard errors
tbl_beh_mod <- tbl_regression(model_beh_mod, exponentiate = TRUE, label = list( 
      ace_counts ~ "ACEs",
      CFI_dich ~ "Child Flourishing Index",
      fpl_category ~ "Federal Poverty Level",
      caregiver_education ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex",
      survey_year ~ "Survey Year",
      caregiver_mental_health ~ "Caregiver Mental Health",
      child_gen_health ~ "Child General Health",
      caregiver_composition ~ "Caregiver Composition",
      physical_resources ~ "Neighborhood Physical Resources",
      environmental_threats ~ "Neighborhood Environmental Threats",
      collective_efficacy ~ "Neighborhood Collective Efficacy",
      special_healthcare_needs ~ "Special Healthcare Needs",
      total_severity_score ~ "NDD Severity (Total)",
      HAI_category ~ "Healthcare Access Index",
      barriers_cat ~ "Access Barriers Index",
      insurance_type ~ "Insurance Type",
      IAI_category ~ "Insurance Adequacy Index"
    )) %>%
  add_q() %>% # adjusts p-values for multiple testing (FDR)
  bold_p(t = 0.05) %>%
  bold_p(t = 0.05, q = TRUE) %>% # now bold q-values under the threshold of 0.05
  italicize_levels() %>%
  gtsummary::as_gt() %>%
  gt::tab_source_note(gt::md("*Note*. This logistic regression model examines the moderating effect of the Child Flourishing Index (CFI) on the association between Adverse Childhood Experiences (ACEs) and behavioral problem outcomes in neurodivergent children. Survey weights from the National Survey of Children's Health (NSCH) were applied to ensure the representativeness of the target population. Robust standard errors (HC0) were used to account for potential heteroscedasticity not fully addressed by the survey weights. HC0 provides the original White standard errors, which are suitable for large samples and do not include small-sample bias corrections. Covariates include child demographics, family characteristics, neighborhood conditions and survey year, with FDR-adjusted p-values to account for multiple comparisons."))

# Print table to review results
tbl_beh_mod

    studentized Breusch-Pagan test

data:  model_beh_mod
BP = 1219.8, df = 52, p-value < 2.2e-16

Characteristic	OR1	95% CI1	p-value	q-value2
ACEs
0 ACEs	—	—
1 ACE	1.04	0.77, 1.41	0.8	>0.9
2-3 ACEs	0.84	0.61, 1.18	0.3	0.8
4+ ACEs	1.19	0.82, 1.73	0.4	0.8
Child Flourishing Index
Not Flourishing	—	—
Flourishing	0.79	0.36, 1.72	0.5	0.9
Child Race
Non-Hispanic White	—	—
Non-Hispanic Black or African American	1.06	0.75, 1.48	0.8	0.9
Non-Hispanic American Indian or Alaska Native	0.80	0.33, 1.92	0.6	0.9
Non-Hispanic Asian, Native Hawaiian, or Pacific Islander	0.61	0.31, 1.21	0.2	0.5
Non-Hispanic Other or Mixed Race	0.93	0.67, 1.29	0.7	0.9
Hispanic/Latino of Any Race	0.90	0.67, 1.22	0.5	0.9
Child Sex
Male	—	—
Female	1.04	0.83, 1.30	0.8	0.9
Federal Poverty Level
<100%	—	—
100%-199%	0.92	0.63, 1.32	0.6	0.9
200%-299%	1.12	0.74, 1.69	0.6	0.9
300%-399%	0.81	0.52, 1.27	0.4	0.8
≥400%	0.90	0.58, 1.39	0.6	0.9
Caregiver Education
Less Than High School	—	—
High School/Vocational	0.83	0.42, 1.63	0.6	0.9
Some College/Associate Degree	0.77	0.39, 1.50	0.4	0.8
College Degree/Higher	0.67	0.34, 1.32	0.2	0.7
Child Age	0.83	0.80, 0.86	<0.001	<0.001
Survey Year
2016	—	—
2017	0.80	0.53, 1.22	0.3	0.8
2018	0.78	0.54, 1.13	0.2	0.6
2019	0.92	0.61, 1.41	0.7	0.9
2020	1.40	0.91, 2.14	0.12	0.5
2021	1.02	0.70, 1.48	>0.9	>0.9
2022	1.07	0.75, 1.52	0.7	0.9
2023	0.74	0.52, 1.06	0.11	0.5
Caregiver Mental Health
Excellent	—	—
Very Good	1.08	0.84, 1.39	0.5	0.9
Good	1.00	0.74, 1.35	>0.9	>0.9
Fair	1.12	0.69, 1.80	0.6	0.9
Poor	1.15	0.34, 3.92	0.8	>0.9
Child General Health
Excellent	—	—
Very Good	1.33	1.05, 1.67	0.018	0.12
Good	1.29	0.96, 1.74	0.10	0.5
Fair	2.49	1.40, 4.40	0.002	0.015
Poor	26.7	5.54, 129	<0.001	<0.001
Caregiver Composition
Dual-caregiver household	—	—
Single-caregiver household	0.90	0.70, 1.14	0.4	0.8
Neighborhood Physical Resources	1.00	0.90, 1.12	>0.9	>0.9
Neighborhood Environmental Threats	0.96	0.85, 1.09	0.5	0.9
Neighborhood Collective Efficacy	1.11	0.99, 1.25	0.084	0.4
Special Healthcare Needs
Non-SHCN	—	—
SHCN	3.23	2.46, 4.23	<0.001	<0.001
NDD Severity (Total)	1.17	1.12, 1.22	<0.001	<0.001
Healthcare Access Index
High Access	—	—
Medium Access	0.67	0.53, 0.84	<0.001	0.007
Low Access	1.03	0.72, 1.48	0.9	>0.9
Access Barriers Index
No Barriers	—	—
Few Barriers	0.94	0.66, 1.33	0.7	0.9
Moderate Barriers	0.74	0.37, 1.51	0.4	0.8
Many Barriers	1.07	0.29, 3.90	>0.9	>0.9
Insurance Type
Private Only	—	—
Private + Public	1.38	0.89, 2.15	0.2	0.5
Public Only	1.07	0.81, 1.41	0.6	0.9
Uninsured	9.45	0.94, 95.0	0.057	0.3
Insurance Adequacy Index
High Adequacy	—	—
Moderate Adequacy	0.70	0.55, 0.90	0.005	0.037
Low Adequacy	0.86	0.62, 1.18	0.4	0.8
ACEs * Child Flourishing Index
1 ACE * Flourishing	0.50	0.11, 2.20	0.4	0.8
2-3 ACEs * Flourishing	2.37	0.61, 9.22	0.2	0.7
4+ ACEs * Flourishing	0.51	0.13, 1.99	0.3	0.8
Note. This logistic regression model examines the moderating effect of the Child Flourishing Index (CFI) on the association between Adverse Childhood Experiences (ACEs) and behavioral problem outcomes in neurodivergent children. Survey weights from the National Survey of Children’s Health (NSCH) were applied to ensure the representativeness of the target population. Robust standard errors (HC0) were used to account for potential heteroscedasticity not fully addressed by the survey weights. HC0 provides the original White standard errors, which are suitable for large samples and do not include small-sample bias corrections. Covariates include child demographics, family characteristics, neighborhood conditions and survey year, with FDR-adjusted p-values to account for multiple comparisons.
1 OR = Odds Ratio, CI = Confidence Interval
2 False discovery rate correction for multiple testing

1. The anxiety model did not yield significant interactions, so we do not proceed with stratified post-hoc analyses. Instead, we will run main effect models without interaction terms.
2. The depression model did not yield significant interactions, so we do not proceed with stratified post-hoc analyses. Instead, we will run main effect models without interaction terms.
3. The behavioral problems model did not yield significant interactions, so we do not proceed stratified with post-hoc analyses. Instead, we will run main effect models without interaction terms.

Main Effect Models

Fit three logistic regression models with survey weights to examine the main effects of ACE counts and the Child Flourishing Index (CFI_dich) on psychopathology outcomes (anxiety, depression, behavioral problems). Extract residuals, check heteroscedasticity using Breusch-Pagan tests, calculate Cook’s distance for influential observations, and obtain robust standard errors. Create summary tables with robust standard errors and False Discovery Rate (FDR) adjusted p-values for each model.

Anxiety Model

# Logistic regression model for main effects using robust standard errors
model_anx <- svyglm(Anxiety ~ ace_counts + CFI_dich + sc_race_cat + sc_sex_cat + 
                        fpl_category + caregiver_education + sc_age_years + survey_year + 
                        caregiver_mental_health + child_gen_health + 
                        caregiver_composition + physical_resources + 
                        environmental_threats + collective_efficacy +
                        special_healthcare_needs + total_severity_score + 
                        HAI_category + barriers_cat + insurance_type + IAI_category, 
                        design = design, family = quasibinomial())

# Extract residuals from the logistic regression model
residuals_anx <- residuals(model_anx, type = "deviance")

# Plot residuals to visually inspect them
plot(residuals_anx, main = "Residuals Plot for Anxiety Model", ylab = "Residuals", xlab = "Index", pch = 20, col = "blue")
abline(h = 0, col = "red", lty = 2)

# Perform the Breusch-Pagan test for heteroscedasticity
bptest(model_anx)

    studentized Breusch-Pagan test

data:  model_anx
BP = 849.67, df = 49, p-value < 2.2e-16

# Calculate Cook's distance for influential observations (leverage) in each model
cooks_anx <- cooks.distance(model_anx)

# Plot for Anxiety model
plot(cooks_anx, type = "p", main = "Cook's Distance for Anxiety Model",
     ylab = "Cook's Distance", xlab = "Observation Index")
abline(h = 4 / length(cooks_anx), col = "red", lty = 2) # Add a reference line for influential points

# Obtain robust standard errors for Anxiety model
robust_anx <- coeftest(model_anx, vcov = vcovHC(model_anx, type = "HC0"))

# Create a summary table for the Anxiety model with robust standard errors
tbl_anx <- tbl_regression(model_anx, exponentiate = TRUE, label = list(
      ace_counts ~ "ACEs",
      CFI_dich ~ "Child Flourishing Index",
      fpl_category ~ "Federal Poverty Level",
      caregiver_education ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex",
      survey_year ~ "Survey Year",
      caregiver_mental_health ~ "Caregiver Mental Health",
      child_gen_health ~ "Child General Health",
      caregiver_composition ~ "Caregiver Composition",
      physical_resources ~ "Neighborhood Physical Resources",
      environmental_threats ~ "Neighborhood Environmental Threats",
      collective_efficacy ~ "Neighborhood Collective Efficacy",
      special_healthcare_needs ~ "Special Healthcare Needs",
      total_severity_score ~ "NDD Severity (Total)",
      HAI_category ~ "Healthcare Access Index",
      barriers_cat ~ "Access Barriers Index",
      insurance_type ~ "Insurance Type",
      IAI_category ~ "Insurance Adequacy Index"
    )) %>%
  add_q() %>% # adjusts p-values for multiple testing (FDR)
  bold_p(t = 0.05) %>%
  bold_p(t = 0.05, q = TRUE) %>% # now bold q-values under the threshold of 0.05
  italicize_levels() %>%
  gtsummary::as_gt() %>%
  gt::tab_source_note(gt::md("*Note*. This logistic regression model examines the association between Adverse Childhood Experiences (ACEs) and anxiety outcomes in neurodivergent children, including main effects of the Child Flourishing Index (CFI) without interaction terms. Survey weights from the National Survey of Children's Health (NSCH) were applied to ensure the representativeness of the target population. Robust standard errors (HC0) were used to account for potential heteroscedasticity not fully addressed by the survey weights. HC0 provides the original White standard errors, which are suitable for large samples and do not include small-sample bias corrections. Covariates include child demographics, family characteristics, neighborhood conditions, and survey year, with FDR-adjusted p-values to account for multiple comparisons."))

# Print table to review results
tbl_anx

Characteristic	OR1	95% CI1	p-value	q-value2
ACEs
0 ACEs	—	—
1 ACE	1.15	0.77, 1.72	0.5	0.7
2-3 ACEs	1.38	0.89, 2.14	0.2	0.3
4+ ACEs	1.56	0.96, 2.54	0.074	0.2
Child Flourishing Index
Not Flourishing	—	—
Flourishing	1.27	0.62, 2.63	0.5	0.8
Child Race
Non-Hispanic White	—	—
Non-Hispanic Black or African American	0.78	0.45, 1.36	0.4	0.6
Non-Hispanic American Indian or Alaska Native	0.58	0.23, 1.44	0.2	0.4
Non-Hispanic Asian, Native Hawaiian, or Pacific Islander	0.33	0.17, 0.65	0.001	0.010
Non-Hispanic Other or Mixed Race	0.95	0.61, 1.46	0.8	0.9
Hispanic/Latino of Any Race	0.61	0.43, 0.86	0.005	0.021
Child Sex
Male	—	—
Female	0.98	0.73, 1.32	0.9	>0.9
Federal Poverty Level
<100%	—	—
100%-199%	0.59	0.35, 1.00	0.050	0.14
200%-299%	0.77	0.43, 1.37	0.4	0.6
300%-399%	0.51	0.28, 0.93	0.028	0.093
≥400%	0.75	0.42, 1.34	0.3	0.6
Caregiver Education
Less Than High School	—	—
High School/Vocational	3.09	1.41, 6.75	0.005	0.021
Some College/Associate Degree	2.97	1.43, 6.16	0.003	0.021
College Degree/Higher	4.20	1.98, 8.92	<0.001	0.003
Child Age	0.99	0.95, 1.04	0.7	0.9
Survey Year
2016	—	—
2017	1.12	0.60, 2.09	0.7	0.9
2018	0.92	0.56, 1.52	0.7	0.9
2019	0.64	0.39, 1.06	0.084	0.2
2020	1.20	0.67, 2.14	0.5	0.8
2021	1.90	1.22, 2.96	0.005	0.021
2022	1.44	0.89, 2.33	0.14	0.3
2023	1.23	0.76, 2.01	0.4	0.6
Caregiver Mental Health
Excellent	—	—
Very Good	1.24	0.90, 1.72	0.2	0.4
Good	1.76	1.21, 2.58	0.003	0.021
Fair	2.43	1.42, 4.17	0.001	0.010
Poor	5.80	1.72, 19.6	0.005	0.021
Child General Health
Excellent	—	—
Very Good	1.04	0.77, 1.39	0.8	0.9
Good	1.59	1.05, 2.40	0.027	0.093
Fair	2.02	0.86, 4.76	0.11	0.2
Poor	1.52	0.39, 5.88	0.5	0.8
Caregiver Composition
Dual-caregiver household	—	—
Single-caregiver household	1.00	0.70, 1.43	>0.9	>0.9
Neighborhood Physical Resources	0.81	0.69, 0.95	0.011	0.042
Neighborhood Environmental Threats	1.01	0.85, 1.20	>0.9	>0.9
Neighborhood Collective Efficacy	0.87	0.74, 1.02	0.092	0.2
Special Healthcare Needs
Non-SHCN	—	—
SHCN	1.78	1.26, 2.53	0.001	0.010
NDD Severity (Total)	1.14	1.08, 1.21	<0.001	<0.001
Healthcare Access Index
High Access	—	—
Medium Access	0.36	0.25, 0.51	<0.001	<0.001
Low Access	0.92	0.53, 1.60	0.8	0.9
Access Barriers Index
No Barriers	—	—
Few Barriers	0.88	0.56, 1.37	0.6	0.8
Moderate Barriers	0.90	0.38, 2.08	0.8	0.9
Many Barriers	0.76	0.20, 2.97	0.7	0.9
Insurance Type
Private Only	—	—
Private + Public	1.06	0.61, 1.84	0.8	0.9
Public Only	0.70	0.46, 1.04	0.078	0.2
Uninsured	0.09	0.01, 0.93	0.044	0.13
Insurance Adequacy Index
High Adequacy	—	—
Moderate Adequacy	1.22	0.85, 1.75	0.3	0.5
Low Adequacy	0.91	0.57, 1.46	0.7	0.9
Note. This logistic regression model examines the association between Adverse Childhood Experiences (ACEs) and anxiety outcomes in neurodivergent children, including main effects of the Child Flourishing Index (CFI) without interaction terms. Survey weights from the National Survey of Children’s Health (NSCH) were applied to ensure the representativeness of the target population. Robust standard errors (HC0) were used to account for potential heteroscedasticity not fully addressed by the survey weights. HC0 provides the original White standard errors, which are suitable for large samples and do not include small-sample bias corrections. Covariates include child demographics, family characteristics, neighborhood conditions, and survey year, with FDR-adjusted p-values to account for multiple comparisons.
1 OR = Odds Ratio, CI = Confidence Interval
2 False discovery rate correction for multiple testing

Depression Model

# Logistic regression model for main effects using robust standard errors
model_dep <- svyglm(Depression ~ ace_counts + CFI_dich + sc_race_cat + sc_sex_cat + 
                        fpl_category + caregiver_education + sc_age_years + survey_year + 
                        caregiver_mental_health + child_gen_health + 
                        caregiver_composition + physical_resources + 
                        environmental_threats + collective_efficacy +
                        special_healthcare_needs + total_severity_score + 
                        HAI_category + barriers_cat + insurance_type + IAI_category,  
                        design = design, family = quasibinomial())

# Extract residuals from the logistic regression model
residuals_dep <- residuals(model_dep, type = "deviance")

# Plot residuals to visually inspect them
plot(residuals_dep, main = "Residuals Plot for Depression Model", ylab = "Residuals", xlab = "Index", pch = 20, col = "blue")
abline(h = 0, col = "red", lty = 2)

# Perform the Breusch-Pagan test for heteroscedasticity
bptest(model_dep)

    studentized Breusch-Pagan test

data:  model_dep
BP = 467.39, df = 49, p-value < 2.2e-16

# Calculate Cook's distance for influential observations (leverage) in each model
cooks_dep <- cooks.distance(model_dep)

# Plot for Depression model
plot(cooks_dep, type = "p", main = "Cook's Distance for Depression Model",
     ylab = "Cook's Distance", xlab = "Observation Index")
abline(h = 4 / length(cooks_dep), col = "red", lty = 2) # Add a reference line for influential points

# Obtain robust standard errors for Depression model
robust_dep <- coeftest(model_dep, vcov = vcovHC(model_dep, type = "HC0"))

# Create a summary table for the Depression model with robust standard errors
tbl_dep <- tbl_regression(model_dep, exponentiate = TRUE, label = list(
      ace_counts ~ "ACEs",
      CFI_dich ~ "Child Flourishing Index",
      fpl_category ~ "Federal Poverty Level",
      caregiver_education ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex",
      survey_year ~ "Survey Year",
      caregiver_mental_health ~ "Caregiver Mental Health",
      child_gen_health ~ "Child General Health",
      caregiver_composition ~ "Caregiver Composition",
      physical_resources ~ "Neighborhood Physical Resources",
      environmental_threats ~ "Neighborhood Environmental Threats",
      collective_efficacy ~ "Neighborhood Collective Efficacy",
      special_healthcare_needs ~ "Special Healthcare Needs",
      total_severity_score ~ "NDD Severity (Total)",
      HAI_category ~ "Healthcare Access Index",
      barriers_cat ~ "Access Barriers Index",
      insurance_type ~ "Insurance Type",
      IAI_category ~ "Insurance Adequacy Index"
    )) %>%
  add_q() %>% # adjusts p-values for multiple testing (FDR)
  bold_p(t = 0.05) %>%
  bold_p(t = 0.05, q = TRUE) %>% # now bold q-values under the threshold of 0.05
  italicize_levels() %>%
  gtsummary::as_gt() %>%
  gt::tab_source_note(gt::md("*Note*. This logistic regression model examines the association between Adverse Childhood Experiences (ACEs) and depression outcomes in neurodivergent children, including main effects of the Child Flourishing Index (CFI) without interaction terms. Survey weights from the National Survey of Children's Health (NSCH) were applied to ensure the representativeness of the target population. Robust standard errors (HC0) were used to account for potential heteroscedasticity not fully addressed by the survey weights. HC0 provides the original White standard errors, which are suitable for large samples and do not include small-sample bias corrections. Covariates include child demographics, family characteristics, neighborhood conditions, and survey year, with FDR-adjusted p-values to account for multiple comparisons."))

# Print table to review results
tbl_dep

Characteristic	OR1	95% CI1	p-value	q-value2
ACEs
0 ACEs	—	—
1 ACE	1.50	1.02, 2.23	0.042	0.13
2-3 ACEs	1.85	1.23, 2.79	0.003	0.027
4+ ACEs	2.17	1.41, 3.33	<0.001	0.006
Child Flourishing Index
Not Flourishing	—	—
Flourishing	1.84	0.96, 3.52	0.065	0.2
Child Race
Non-Hispanic White	—	—
Non-Hispanic Black or African American	0.76	0.46, 1.24	0.3	0.5
Non-Hispanic American Indian or Alaska Native	0.77	0.26, 2.29	0.6	0.8
Non-Hispanic Asian, Native Hawaiian, or Pacific Islander	1.29	0.63, 2.65	0.5	0.7
Non-Hispanic Other or Mixed Race	1.45	0.95, 2.22	0.085	0.2
Hispanic/Latino of Any Race	0.58	0.39, 0.88	0.010	0.056
Child Sex
Male	—	—
Female	1.42	1.06, 1.91	0.019	0.083
Federal Poverty Level
<100%	—	—
100%-199%	1.17	0.71, 1.91	0.5	0.7
200%-299%	0.79	0.43, 1.46	0.5	0.7
300%-399%	0.85	0.42, 1.73	0.7	0.8
≥400%	0.98	0.52, 1.84	>0.9	>0.9
Caregiver Education
Less Than High School	—	—
High School/Vocational	1.15	0.55, 2.44	0.7	0.8
Some College/Associate Degree	0.89	0.44, 1.81	0.7	0.8
College Degree/Higher	0.93	0.42, 2.07	0.9	0.9
Child Age	1.00	0.95, 1.05	>0.9	>0.9
Survey Year
2016	—	—
2017	1.46	0.71, 3.00	0.3	0.5
2018	1.91	1.04, 3.51	0.037	0.13
2019	1.14	0.63, 2.07	0.7	0.8
2020	1.94	1.02, 3.67	0.042	0.13
2021	1.49	0.80, 2.78	0.2	0.5
2022	1.45	0.82, 2.56	0.2	0.5
2023	1.35	0.77, 2.38	0.3	0.5
Caregiver Mental Health
Excellent	—	—
Very Good	1.05	0.74, 1.48	0.8	0.8
Good	1.95	1.31, 2.90	0.001	0.010
Fair	2.08	1.09, 3.98	0.027	0.10
Poor	3.89	1.23, 12.3	0.021	0.085
Child General Health
Excellent	—	—
Very Good	0.92	0.68, 1.25	0.6	0.8
Good	1.71	1.14, 2.56	0.010	0.056
Fair	1.57	0.78, 3.16	0.2	0.5
Poor	8.48	1.99, 36.2	0.004	0.027
Caregiver Composition
Dual-caregiver household	—	—
Single-caregiver household	0.59	0.43, 0.81	<0.001	0.010
Neighborhood Physical Resources	0.93	0.78, 1.12	0.5	0.7
Neighborhood Environmental Threats	0.94	0.76, 1.16	0.6	0.7
Neighborhood Collective Efficacy	0.97	0.82, 1.15	0.7	0.8
Special Healthcare Needs
Non-SHCN	—	—
SHCN	2.14	1.47, 3.12	<0.001	0.002
NDD Severity (Total)	1.08	1.02, 1.16	0.016	0.078
Healthcare Access Index
High Access	—	—
Medium Access	0.32	0.22, 0.46	<0.001	<0.001
Low Access	0.73	0.47, 1.12	0.2	0.4
Access Barriers Index
No Barriers	—	—
Few Barriers	0.84	0.60, 1.18	0.3	0.5
Moderate Barriers	0.78	0.34, 1.76	0.5	0.7
Many Barriers	0.78	0.21, 2.93	0.7	0.8
Insurance Type
Private Only	—	—
Private + Public	1.23	0.70, 2.17	0.5	0.7
Public Only	1.07	0.71, 1.61	0.7	0.8
Uninsured	0.23	0.03, 1.58	0.13	0.3
Insurance Adequacy Index
High Adequacy	—	—
Moderate Adequacy	1.17	0.87, 1.58	0.3	0.5
Low Adequacy	1.25	0.76, 2.04	0.4	0.6
Note. This logistic regression model examines the association between Adverse Childhood Experiences (ACEs) and depression outcomes in neurodivergent children, including main effects of the Child Flourishing Index (CFI) without interaction terms. Survey weights from the National Survey of Children’s Health (NSCH) were applied to ensure the representativeness of the target population. Robust standard errors (HC0) were used to account for potential heteroscedasticity not fully addressed by the survey weights. HC0 provides the original White standard errors, which are suitable for large samples and do not include small-sample bias corrections. Covariates include child demographics, family characteristics, neighborhood conditions, and survey year, with FDR-adjusted p-values to account for multiple comparisons.
1 OR = Odds Ratio, CI = Confidence Interval
2 False discovery rate correction for multiple testing

Behavioral Problems Model

# Logistic regression model for main effects using robust standard errors
model_beh <- svyglm(Behavioral_Problems ~ ace_counts + CFI_dich + sc_race_cat + 
                        sc_sex_cat + fpl_category + caregiver_education + sc_age_years + 
                        survey_year + caregiver_mental_health + child_gen_health + 
                        caregiver_composition + physical_resources + environmental_threats +
                        collective_efficacy + special_healthcare_needs + 
                        total_severity_score + HAI_category + barriers_cat + 
                        insurance_type + IAI_category,
                        design = design, family = quasibinomial())

# Extract residuals from the logistic regression model
residuals_beh <- residuals(model_beh, type = "deviance")

# Plot residuals to visually inspect them
plot(residuals_beh, main = "Residuals Plot for Behavioral Problems Model", ylab = "Residuals", xlab = "Index", pch = 20, col = "blue")
abline(h = 0, col = "red", lty = 2)

# Perform the Breusch-Pagan test for heteroscedasticity
bptest(model_beh)

    studentized Breusch-Pagan test

data:  model_beh
BP = 1228.6, df = 49, p-value < 2.2e-16

# Calculate Cook's distance for influential observations (leverage) in each model
cooks_beh <- cooks.distance(model_beh)

# Plot for Behavioral Problems model
plot(cooks_beh, type = "p", main = "Cook's Distance for Behavioral Problems Model",
     ylab = "Cook's Distance", xlab = "Observation Index")
abline(h = 4 / length(cooks_beh), col = "red", lty = 2) # Add a reference line for influential points

# Obtain robust standard errors for Behavioral Problems model
robust_beh <- coeftest(model_beh, vcov = vcovHC(model_beh, type = "HC0"))

# Create a summary table for the Behavioral Problems model with robust standard errors
tbl_beh <- tbl_regression(model_beh, exponentiate = TRUE, label = list( 
      ace_counts ~ "ACEs",
      CFI_dich ~ "Child Flourishing Index",
      fpl_category ~ "Federal Poverty Level",
      caregiver_education ~ "Caregiver Education",
      sc_age_years ~ "Child Age",
      sc_race_cat ~ "Child Race",
      sc_sex_cat ~ "Child Sex",
      survey_year ~ "Survey Year",
      caregiver_mental_health ~ "Caregiver Mental Health",
      child_gen_health ~ "Child General Health",
      caregiver_composition ~ "Caregiver Composition",
      physical_resources ~ "Neighborhood Physical Resources",
      environmental_threats ~ "Neighborhood Environmental Threats",
      collective_efficacy ~ "Neighborhood Collective Efficacy",
      special_healthcare_needs ~ "Special Healthcare Needs",
      total_severity_score ~ "NDD Severity (Total)",
      HAI_category ~ "Healthcare Access Index",
      barriers_cat ~ "Access Barriers Index",
      insurance_type ~ "Insurance Type",
      IAI_category ~ "Insurance Adequacy Index"
    )) %>%
  add_q() %>% # adjusts p-values for multiple testing (FDR)
  bold_p(t = 0.05) %>%
  bold_p(t = 0.05, q = TRUE) %>% # now bold q-values under the threshold of 0.05
  italicize_levels() %>%
  gtsummary::as_gt() %>%
  gt::tab_source_note(gt::md("*Note*. This logistic regression model examines the association between Adverse Childhood Experiences (ACEs) and behavioral problem outcomes in neurodivergent children, including main effects of the Child Flourishing Index (CFI) without interaction terms. Survey weights from the National Survey of Children's Health (NSCH) were applied to ensure the representativeness of the target population. Robust standard errors (HC0) were used to account for potential heteroscedasticity not fully addressed by the survey weights. HC0 provides the original White standard errors, which are suitable for large samples and do not include small-sample bias corrections. Covariates include child demographics, family characteristics, neighborhood conditions, and survey year, with FDR-adjusted p-values to account for multiple comparisons."))

# Print table to review results
tbl_beh

Characteristic	OR1	95% CI1	p-value	q-value2
ACEs
0 ACEs	—	—
1 ACE	1.00	0.74, 1.37	>0.9	>0.9
2-3 ACEs	0.86	0.62, 1.19	0.4	0.8
4+ ACEs	1.18	0.81, 1.71	0.4	0.8
Child Flourishing Index
Not Flourishing	—	—
Flourishing	0.66	0.27, 1.59	0.4	0.8
Child Race
Non-Hispanic White	—	—
Non-Hispanic Black or African American	1.06	0.76, 1.49	0.7	0.9
Non-Hispanic American Indian or Alaska Native	0.81	0.33, 1.95	0.6	0.9
Non-Hispanic Asian, Native Hawaiian, or Pacific Islander	0.61	0.31, 1.20	0.2	0.5
Non-Hispanic Other or Mixed Race	0.95	0.68, 1.33	0.8	0.9
Hispanic/Latino of Any Race	0.90	0.66, 1.22	0.5	0.9
Child Sex
Male	—	—
Female	1.04	0.83, 1.31	0.7	0.9
Federal Poverty Level
<100%	—	—
100%-199%	0.91	0.63, 1.31	0.6	0.9
200%-299%	1.11	0.74, 1.68	0.6	0.9
300%-399%	0.80	0.51, 1.25	0.3	0.8
≥400%	0.88	0.57, 1.36	0.6	0.9
Caregiver Education
Less Than High School	—	—
High School/Vocational	0.83	0.43, 1.63	0.6	0.9
Some College/Associate Degree	0.77	0.40, 1.51	0.5	0.9
College Degree/Higher	0.68	0.34, 1.35	0.3	0.8
Child Age	0.83	0.80, 0.86	<0.001	<0.001
Survey Year
2016	—	—
2017	0.81	0.53, 1.24	0.3	0.8
2018	0.79	0.54, 1.14	0.2	0.6
2019	0.92	0.60, 1.40	0.7	0.9
2020	1.38	0.91, 2.11	0.13	0.5
2021	1.01	0.70, 1.47	>0.9	>0.9
2022	1.06	0.74, 1.50	0.8	0.9
2023	0.74	0.52, 1.06	0.10	0.5
Caregiver Mental Health
Excellent	—	—
Very Good	1.08	0.83, 1.40	0.6	0.9
Good	0.99	0.73, 1.34	>0.9	>0.9
Fair	1.11	0.69, 1.79	0.7	0.9
Poor	1.14	0.33, 3.88	0.8	>0.9
Child General Health
Excellent	—	—
Very Good	1.33	1.05, 1.69	0.016	0.10
Good	1.30	0.96, 1.75	0.090	0.4
Fair	2.51	1.41, 4.44	0.002	0.013
Poor	26.7	5.54, 128	<0.001	<0.001
Caregiver Composition
Dual-caregiver household	—	—
Single-caregiver household	0.90	0.70, 1.15	0.4	0.8
Neighborhood Physical Resources	1.00	0.90, 1.11	>0.9	>0.9
Neighborhood Environmental Threats	0.96	0.85, 1.09	0.5	0.9
Neighborhood Collective Efficacy	1.10	0.98, 1.24	0.12	0.5
Special Healthcare Needs
Non-SHCN	—	—
SHCN	3.20	2.44, 4.20	<0.001	<0.001
NDD Severity (Total)	1.17	1.12, 1.22	<0.001	<0.001
Healthcare Access Index
High Access	—	—
Medium Access	0.68	0.54, 0.85	<0.001	0.009
Low Access	1.04	0.73, 1.49	0.8	>0.9
Access Barriers Index
No Barriers	—	—
Few Barriers	0.94	0.66, 1.33	0.7	0.9
Moderate Barriers	0.74	0.37, 1.50	0.4	0.8
Many Barriers	1.06	0.29, 3.88	>0.9	>0.9
Insurance Type
Private Only	—	—
Private + Public	1.38	0.89, 2.16	0.2	0.5
Public Only	1.06	0.81, 1.40	0.7	0.9
Uninsured	9.49	0.95, 95.2	0.056	0.3
Insurance Adequacy Index
High Adequacy	—	—
Moderate Adequacy	0.70	0.54, 0.89	0.005	0.032
Low Adequacy	0.85	0.62, 1.17	0.3	0.8
Note. This logistic regression model examines the association between Adverse Childhood Experiences (ACEs) and behavioral problem outcomes in neurodivergent children, including main effects of the Child Flourishing Index (CFI) without interaction terms. Survey weights from the National Survey of Children’s Health (NSCH) were applied to ensure the representativeness of the target population. Robust standard errors (HC0) were used to account for potential heteroscedasticity not fully addressed by the survey weights. HC0 provides the original White standard errors, which are suitable for large samples and do not include small-sample bias corrections. Covariates include child demographics, family characteristics, neighborhood conditions, and survey year, with FDR-adjusted p-values to account for multiple comparisons.
1 OR = Odds Ratio, CI = Confidence Interval
2 False discovery rate correction for multiple testing

Session Information

---

To enhance reproducibility, details about the working environment used for these analyses can be found below.

sessionInfo()

R version 4.4.1 (2024-06-14)
Platform: aarch64-apple-darwin20
Running under: macOS 15.1.1

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.0

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Chicago
tzcode source: internal

attached base packages:
[1] grid      stats     graphics  grDevices utils     datasets  methods  
[8] base     

other attached packages:
 [1] ggdist_3.3.2     lubridate_1.9.3  forcats_1.0.0    stringr_1.5.1   
 [5] purrr_1.0.2      tibble_3.2.1     tidyverse_2.0.0  ggpubr_0.6.0    
 [9] crosstable_0.8.1 psych_2.4.6.26   gt_0.11.1        ggeffects_1.7.0 
[13] gtsummary_2.0.3  sandwich_3.1-0   lmtest_0.9-40    zoo_1.8-12      
[17] survey_4.4-2     survival_3.6-4   Matrix_1.7-1     sjmisc_2.8.10   
[21] sjlabelled_1.2.0 ggplot2_3.5.1    kableExtra_1.4.0 knitr_1.48      
[25] readr_2.1.5      tidyr_1.3.1      dplyr_1.1.4      pacman_0.5.1    

loaded via a namespace (and not attached):
 [1] DBI_1.2.3               mnormt_2.1.1            gridExtra_2.3          
 [4] rlang_1.1.4             magrittr_2.0.3          compiler_4.4.1         
 [7] systemfonts_1.1.0       vctrs_0.6.5             quadprog_1.5-8         
[10] pkgconfig_2.0.3         crayon_1.5.3            fastmap_1.2.0          
[13] backports_1.5.0         labeling_0.4.3          utf8_1.2.4             
[16] PupillometryR_0.0.5     rmarkdown_2.27          markdown_1.13          
[19] tzdb_0.4.0              haven_2.5.4             ragg_1.3.2             
[22] bit_4.0.5               xfun_0.45               labelled_2.13.0        
[25] jsonlite_1.8.8          uuid_1.2-0              broom_1.0.6            
[28] parallel_4.4.1          R6_2.5.1                stringi_1.8.4          
[31] car_3.1-2               Rcpp_1.0.13-1           base64enc_0.1-3        
[34] splines_4.4.1           timechange_0.3.0        tidyselect_1.2.1       
[37] rstudioapi_0.16.0       abind_1.4-5             yaml_2.3.9             
[40] lattice_0.22-6          withr_3.0.2             flextable_0.9.7        
[43] askpass_1.2.0           evaluate_0.24.0         zip_2.3.1              
[46] xml2_1.3.6              pillar_1.9.0            carData_3.0-5          
[49] checkmate_2.3.1         insight_0.20.5          distributional_0.5.0   
[52] generics_0.1.3          vroom_1.6.5             hms_1.1.3              
[55] commonmark_1.9.1        munsell_0.5.1           scales_1.3.0           
[58] glue_1.8.0              gdtools_0.4.0           tools_4.4.1            
[61] data.table_1.15.4       ggsignif_0.6.4          mitools_2.4            
[64] cards_0.3.0             colorspace_2.1-0        nlme_3.1-164           
[67] cli_3.6.3               textshaping_0.4.0       officer_0.6.7          
[70] fontBitstreamVera_0.1.1 fansi_1.0.6             broom.helpers_1.17.0   
[73] viridisLite_0.4.2       svglite_2.1.3           gtable_0.3.5           
[76] rstatix_0.7.2           sass_0.4.9              digest_0.6.36          
[79] fontquiver_0.2.1        htmlwidgets_1.6.4       farver_2.1.2           
[82] htmltools_0.5.8.1       lifecycle_1.0.4         fontLiberation_0.1.0   
[85] openssl_2.2.0           bit64_4.0.5